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Related Application 

This application claims the benefit of U.S. Provisional Patent Application No.. 
60/370,507, filed April 4, 2002, entitled "3D Digital Library System," which is incorporated 
herein by reference. 

Background 

This invention relates generally to information storage and retrieval. More specifically, 
it relates to a system and method for the storing, archiving, query and analysis of information 
relatiiig to three-dimensional (3D) materials and objects. The system can be used within a 
distributed environment to catalog, organize and support interaction, wifii 3I> material and 
objects at a feature-based level. 

The understanding of 3D structures is essential to many disciplines, including scientific 
disciplines sudh as archaeology and physical anthropology. For example, (1) arefcaeologists 
study the 3D form of Native American pottery to characterize the development of prehistoric 
cultures; (2) lithic research seeks to understand the development of stone tool technology over 
time by conjoining the remams of individual lithic reduction sequences; (3) comparative 
quantitative analysis of the 3D joint sutfeces of human and nonhumaa prinaate bones and the 
effects of surface area, cxirvature, andcongruency on the mechaiiics;.of a joint can fiiither the 
understanding of physical anthropologists who focus xm the Eeconstroettpn of manipuladve and 
locomotor behavior of early humans. To date, the nietho<J§ that aothrppolpgists use to measure, 
record and classi^if 3D data s^uch, as., c^d^peis, ppolx^ss^h^ Ctfas^iiigSji. and ;.Sie use of 
classiiScation schemes that reduce dbjecte fo,siai|>Ii^e4-4^^::^ta!ES ^ inaicteipiate suxwsiitiily 
represraiti^g such complex data, • . . 

Within the past decade, technologies for c^turing SD objects hssie i^ex^ed and 
iinprbved dramatically^ offering the pos^^ obtain and stgr^accikafen^^ 
of objects in digital form. Notwithstand^; tiias jptojgE^ 3D 
data representation, storage, retrieval, semaritie amchineriti ^i^^fi^ 
eSlraction and developmert of a yisM queiy s3rat€to^^ 
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Summary 

Jd. accordance with the invention, there is provided a computer syst^ add method for 
the storage, archiving, query and retrieval of infonnation relating to 3I>objects. The system 
includes data acquisition means for acquiring point coordinate data about a three-dimensional 
object, a database component, a processor and a user interface. The processor is operable to 
generate modeled data fix>m the point coordmate data and to segment the modeled data into 
feature data represaatmg a plurality of features of the object The data is organized so that 
features of the 3D objects can be automatically extracted for online query and retrieval. The 
processor is opCTable to store the modeled data and the feature data in the database component 
and retrieve modeled data and feature data from the database component using search criteria 
representing object "features of interest. The user interface is operative with the processor to 
allow a user input search criteria to processor and to display data retrieved by the processor as a 
representation of an object feature. 

The point coordinate data can be surface data or volume data. The modeled data can 
comprise 31) triangular mesh data segmented into subsets of vertices. Extracted features can 
include points, curves, surface data and volume data. Segmentation of the modeled data is 
-perfoiemed to simplify the data and to organize the data into features or regions that are 
TrieaningM for applications or to users. One, advantageous method for segmenting the modeled 
data uses a watershed segmentation method with in^rbved curvature estimation. Another 
•advantageous method uses a hybrid segmentation scheme. Data compression can be used to 
compress the modeled data. An advantageous method for compressing the modeled data uses B- 
spliae curves. Another advantag^otis. method uses subdivision sui£ice con^ression. Two 
advantageous methods for isegmentiiig vblume use WeibuU E-SD fields and Greedy comiected 
xoiopoHeiit M>el&g refitii^Eneal 

Hie user interfece includes a graphic user interface that provide a visual quay system 
ssad method &at allows & skefcEt'tiased seiorah of tfie diabase. Repres^tative object slia^ 
' bb:«dec^-^om a pafis^ aad tbdSSed, of a ^£)nii 'proJSIe sketch can be cre^ed using ^ 
visuM qnery. lie initial qiieiy request can be nuide in a vanety of modes ineliiding text, vector 
gcsp^^ ^ aad 31> models. The user iaierface also can include text and 

numeric fields for p arallel query of descriptive aBd derived data vidthin the' databases. The visual 
que^ interface displays seapch results, which can include textiaal data and i^a^lac data, 
-iaoluding 2D and 3D rq»r^entations of objects that niatch^^t^^ 

Thz systmi md Tii^ ijsed vwtii models fo ajUiress ^various Ijpes of data 

selecfed to r^reseaf diffefeiit design diatleziges, volutnetdc and i^»afeil elemoits, ^aidmodeiiug 
re(]uireaneiit& . 

3 - •' . - ■ -•" • • 
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Brief Description Of The Draivings 

The accompanying drawings, which are incorporated in and consfittEte' part of the 
specification, illustrate the presently preferred embodiments and methods of the invention 
Togdher with the general description given above and flie detailed description of the preferred 
embodiments and methods given below, they serve to explain the principles of the invention. 

FIG. 1 is a functional block diagram of an exemplary computer for system storing, 
archiving, query and retrieval of information relating to 3D objects in accordance wifli tlie 
invention. 

FIG. 2 shows an example of an XiSlL schana according to the present invention. 

FIG. 3 shows an exemplary region editor screen- display for editing modeled data 
representing 3D objects. 

FIG. 4 shows an exemplary query input screen display for iiiputting sketch-based, 
numfflic, and text-based search criteria. 

FIG. 4 is a flovvchart of the visual queiy process. 

FIG. 6 shows an exemplary screen for displaying results of a query in the form of 
tliumbnail images and summary data for 3D objects returned by the search. , ■ 

FIG. 7 illustrates an exemplary screen for displajnng further detail of a specific 3D -object 
"selected from among the search results of FIG. 6. 

FIG. 8 is a diagFam illustrating data flow, knowledge absfractibn and qiiery according to 
the invention. 

FIG. 9 illustrates an example of a triangulated, highly decimated point cloud data set 
representing a surface of an object. 

FIG. 10 shows in more detail the data flow and stmcture of the geometric modeling 
module of file systemofHG. I . ' 

FIG. 11 -shows in more detail the data ftow and stracfure of the feature cactractiott, . 
soMya^^md index geseradoa saodxde^nIfhBi system df HO. 

FIG.J21iasbeiSi'lfit^onattjfbiaitte^ , .. . > . 

WG. 13t shows #svo pieces of an oti jec^ Le. a coFe^add'ffafef Hf ai stone ^rtifed, which 
pieces, fit. together and eadi of which pieces has teea segm^ated to ejctract its sdifiuse features 
for comparison. . 

FIG. 14 ishows ciurvature!>asedsqgmentatioil 

FIG, 16 :^ows a <x»apmson of (A) a Gaossiai curvature 
plot ■ ■ ■ ■ ■ 

FIG. 17 shows (A) the nmi-tim box for afbitraJry axes aiid <B) (he min-m^ box for the . 
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axes of the norm ellipse of the points. 

Figure 18: Figure showing local surface fit for curvature estiniafitto 

Figure 19: Considering the curvature sign may yield an incorrect segmentation. 

FIG. 20 shows an exaii^)le where a triangle T,- will be such a gray area: (A) shows an 

example of segm^entation using the watershed method with no hard boundaries; and (B) shows 

the segmentation example of FIG. 20A using the watershed metiiod with such a boundary 

solution. 

FIG. 21 illustrates the creation of triangles on boundaries of an original mesh, according 
to the hybrid method. ♦ 

FIG. 22 shows a mesh segmented using the watershed method and illustrating no hard 
boundary. 

FIG. 23 shows feature vertices and edges of a mesh. 

FIG. 24 shows the steps of t he hj^rid segmentation method. . 

FIG. 25 shows creation of triangles on feature edges of an original mesh. 
FIG- 26 shows assigning regions for the case having multiple labels with iiiultiple 
common labels. 

FIG. 27 shows segmentation of mechanical parts using the hj'brid method. 
FIG. 28 shows segmentation of a turbine and a wheel using tlie hybrid method 
FIG. 29 shows segmentation of a bone using the hybrid method. 

Fig. 30 shows mesh traversal and parallelogram prediction method for Edgebreaker 
algorithm. 

Fig. 31 shows mMesh traversal and recording of ttie error vector in the B -Spline 
algorithiQ. 

Fig, 32 shows Left: nqimai case for recreating the jieometr^ fflid topology. Right: 
haiidling of special case if tbe vertex has beea-vidtedpT.iecx)]^ 

Pig. 33 shows: Exsmple oi O&w. Left:. T^i^le, ]^fesh .ispresentatkoL Bight:, B'^SpIme 
j^pioximadon to the edge iiudpoii]ls.JEac^ .. 

Fig, 34 shows: StaBfoid.Bnaay Top left: .'Fmngle Mesh rq|)PMatitatiom Top ri^: 
Zooiiied in view of the mesh-Bpttpm Ie|t: B..<S|pliQe ^^jpioximation to fbe^^e mid poiiitsv' 
/Bottoni ng^it Zoomed in view of the B-Splme 
^^^^ ^ H ' 

Fig. 36 shows a remiMhiBg process according to one aspect of the invention. ^ 
FIG. 37 shows a flow chart of (herenieshing process of FIG. 36. 

FIG. 38 shows In the Ix)op st^diviaon, level produces two types of . , 

y&&i^: vertex paints edge p^i^^ 
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WCi.3y shows 

FIG. 40 shows the construction and comparison of a number of wetl=iaiown.data sets. 
FIG. 41 shows (A) Polygonal Mesh with watershed defined areas for complete vessel 
and (B) Polygonal Mesh with watershed defined areas for partial vessel. 

FIG. 42 shows shows feature segmentation of complex sh^ed vessel. 
FIG. 43 shows descriptive elements of ceramic vessel. 

FIG. 44 illustrates four kinds of vessel bases, i.e.: (a) convex base; (b) flat base with zero 
. curvature; (c) concave base; and (d) composite base. 

FIG. 45 shows a schema used for archaeological vessels. ,<.. 
FIG. 46 describes a schema used for Utiiics. 

FIG. 47 shows a distribution experiment of a LCM volume data: (a) A real LCM volume 
data; (b) The distribution of the data 

FIG. 48 shows the Weibull distribution with different shape parameters a and scale 
parameter b=l. 2 and vo=0. 

FIG. 49 shows function /B(t,t)/2. It reaches maximum 0J2 at t=0;42. 1 2 
FIG. 50 illustrates classification of noise in a box by using two paratheteis: si and s2, 
(A) Noraial sample case; {64,52, 64, 46,50,54,62,43}, where si =4.475 and s2 =4.9; (B) High 
noise case {240,52, 64, 46,50,54,62,43}, resulted fiiom case (A) just added one high noise term 
on purpose, where si =0.825 and s2=3.6; (C) Low noise case {6,52, 64, 46,50,54,62,43}, 
resulted from case (a) just added one low noise term on purpose,, where « =3 .45 and «H).575; (D) 
Low noise case {6,52, 64, 46,50,54,62,234}, resulted from ease (A) just added one low noise 
term and one high noise term on purpose, where ,1 =0.725 and c ==0.525. 

FIG. 51 shows randomly generating Weibull distnbution using (7) with the scale 
parameter b=1.2 and three distinct shape parameters 8^0,75, a-3 and a=10. (a) The whole 
volume data (b) the partial volume data with two co-centric spheres 

FIG. 52 shows geaeradng.SD Wdboll distributed volume 4a&i dis&ict s&^e and 
scale paiamet«s. In (A^, the blue part (cube) with a=0.75 and b=lJ?^-ihe.ted part (outside 
sphere) with a=3.0 and b=I6.58, and flie green sphere ^nside) with a=l0.ff Mid.b=63.58. The 
naages are rendered using Ray-Casting. -. , 

FIG. 53 shows E-SD plot of 3D control volume data above; Hie coioepart on the lefl 
corresponds to cube outside of spheres, the disp at the nsiddle to tbie.e:stQnia(t spher^j the small 
region on the ri|Jit-top to the internal sphere.. (i^) Cube .gut of sphjaes <b) external sphere (e) 
infernal sphere. : 

FIG. 54 shows the segmentation of control 3D WaljuU dislrilmted Volume data.^ 
coIorisRGB(155,ViO),v*ejievistiievahieattiif poiM<k,yiZ^^ . .. . a:. 
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FIG. 55 shows test bed for our Weibull E-SD modeling scheme. They come from two 

different experiments. The size of image (A) is 512x512x124, the gray4e\^eI-is-from 0 to 255, 

and there is only one cell, its spindle is at the up part of flie cell. The. size of image (B) is 

512x512x146, has the same gray-level with (A). 

FIG. 56 shows the E-SD plots of image (A) and (B) of Fig 55, respectively. The colors 

•have the same means as FIG. 53. In (A), the size of -window is (178, 237; 4, 27) corresponding 

to the segmentation (A) in FIG. 57, and the threshold e T=3A, so there is a blank on the left side. 

There are two clear fingers. In (B), the size of window is (167, 255; 2, 30) corresponding to flie 

segmentation (B) in Fig 57, and tiie threshold e T=34, so tiiere is a blank on the left side. Tliere 

are two clear fingers. The box size is 2x2x2. 

FIG. 57 shows the spindle segmentation of FIG. 55, respectively. The color is RGB(0, 

v,0), wliere is the value at the point (x, y, z). It are rendered by using OpenGL point clouds. 

The segmentation in (A) includes 3902 boxes, and 1 1308 boxes in (B). 

FIG. 58 shows the region corresponding to the small finger in FIG. 58A. It is rendered by 

using OpenGL point clouds from two different views. 

FIG. 60 shows a test bed for our volumetric segmentation scheme. FIG. 60A shows an 

LS CM image of GFAP-labeled astrocytes and DAPI stained nuclei, in whifch the green 

highlights regions targeted with GFAP labeled astrocytes, and ftie red tegions identify DAPI, 

which consistently targets nudei cells. FIG. 60B shows a LSCM image of GFAP-labeled 

astrocytes and electrode implant, in which the red targets implant and the gre^ GFAP. Boft are 

rendered by MIP in 3D- 

FIG. 61 is a plot of fiinction g(a) defined by equation (3) fbr four samples whidi are 

generated.by We&uU pdf with 5 =i= 1:2, and vO = 0:0 and a = 0:7; 2; 3, or 10, respectively. And 

0.7, 2.1, 3.05 and 9.8 are the roots of g^Ol 0 iCone^cHidnig fo fhie samples. 

KEG. 62 shows the WabuU paramet^ a-6 histogram for the LSCM image shown in 

Figure 1(a) and the color legMi^wJiate N(a;b) is t&9 niaa&er of voxels "wiiicli satisfy WefbuU 

parameters^are and 5 

FIG. 63 shows an example of hierarch3rftame for a 2D image, 

FIG. 64 shows segmentation (a) before GCCL, and (b)afl^ GCCL. Both are rendered by 
point clouds in 3D. 

FLQ. 65 iMiistEates sa3aoofhjng connected componenSs using convolution: a GFAP labeled 
astrocyte (a) b6fore smooth&g, aiMi (b) itfier smoottjing^ imd; a nudd (e) before smoothing aad 
(d) after smoothing 

FIG. 66 iBostiatesthe' segnientation lesates^^^^^ LSGM image of flie GFAP-labeled 
astrocytes ami DAPI stained nudei (FIG. 60A) •' 
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FIG. 67 shows the connection analysis for GFAP labeled astroc3^es in LSCM image 

FIG.60A 

FIG. 68 illustrates the segmentation results of the GFAP and implant LSCM image 
Figure 1(b), (a) two chaimel data, and (b) the GFAP-labeled astrocytes, and (c) the electrode 
implant segmented 

Description 

Reference will now be made in more detail to the presently preferred embodiments and 
methods of the invention as illustrated in the accompanying drawings, in which like numerals 
refer to hke parts throughout flie several views. 

The present invention provides a system and method for the storing, archiving, query and 
analysis of information relating to 3D objects. The system can be used within a distributed 
environment to catalog, organize and support interaction with 3D materials and objects. 

The Computet Nehvork System 
FIG. 1 illustrates in schematic block diagram form an exemplary computer netvs'ork 
system 100 for practicing the invention. The computer network system 100 includes a server 
computer system 102 and a client computer system 104, which are coimected by data 
connections 105, 107 to a computer network 106, e.g., an intranet, the Intemet and/or the World 
Wide Web, so that the cliait computer system 104 and the servCT computea: sj'stem 102 can 
communicate. As will be readily apparent to persons skilled in die art, the client cosjputer 
^tem 104 is intended to be representative of a plurality of client computer systems 104^ each of 
which may communicate with the server computer system. tOl via. the netwoik 106, wheffiier . 
siequentialiy or simxiltaneously. It is to be understood fliat the terms cjient and server..|ire to be 
coDstnied in the bix)adest sense, and that aH such constructions of (he terms are int^epded to fall 
Mffi& the s(»ije of the a^^ 
■ ' "' ^TOtii reference to iiis netwoik 106, it is contenq)lated accordrng to the present invention 
tiitk tsafiftof the cHcMts lb4 may aceess andi cMnmwuicaie wiQt fiie network 106 fittou^;.?py data 
'i^tttttiiilriic^bii technblo^^^ For example, the client 104 may comprise, one or more personal 
xdiBputers that' are pM area network (LAN) and the data connections 

105, 107 may be provided by dedicated data lines. Personal Commimieation: Systems (PCS), 
itncibwave, s telg>hones, or any pthep suifaible means. For example, 

the client cbmjjuter 104 tiiat can be connected to the internet via a phone network, such as those 
provided by a local or re^onal tdephoae operating con^aay, or a coaxial cable line. In any 
case, c&dii{si 104 are ad^ted to coinmunicate with flie netwoik. 10^ such flat ip&iicnation may 
fce^ tiaosmitfetif fo and ikom fee ijerver 102, Cig,, thioaji^ <«ie ions; i:!3iJters,.wide area networks 
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(WANs), satellites, hubs, repeaters, bridges and gateways, as is known in the art. Data 

transmissions are typically passed from network to network in packets tbatinciude not only the • 

substantive a^ects of the data transmission, but addresses, error checkiiig information and the 

like. 

The computer network system 10 advantageously makes use of standard Intemet 
protocols including TCP/IP and HTTP. TCP/TP is a common transport layer protocol used by a 
worldwide network of computers. HTTP is a known plication protocol that provides users 
access to files (which can be in different formats such as text, graphics, images, sound, video, 
etc.) using a standard page description language known as Hypertext Markup Language 
(HTML), DHTML or XML. Known HTML web browsers allow for graphical user interface 
(GUI) based access to HTML documents accessible on servers communicatively linked to the 
client 104. These documents are commonly referred to as "web pages." Although the client 104 
and the server cornputer 102 are coupled together via the Intemet, the invention may also be 
implemented over otlier public or private networks or may be employed through a direct 
connection and any such cotmnunicafion implementation is contemplated as falling within the 
scope of the present invention. 

Server Computer System 
The server computer system 102, which has a conventional architecture, includes: a 
central processing unit (CPU), which may be in:q)lemented with a conventional microprocessor^ 
Means for temporary storage of information, which may be implemented with random access 
memory (RAM)j and means for permanent storage of information, which may be inoplemented 
Wiffi read only memory (ROM); and means for mass storage of information, whiph may be 
itaplenienf ed by hard drive or any o&er suitable means of mass storage known in the art, 

If vrfll bb obvious to someone of ordinaiy skill in flie art flia^ 
a variety of other system architectures. As described herein, the exeiiiplary s^tem architecture 
is'for descriptive purposes bi^y. Altlbough the desedpticm ma^ ve^ ^& Um^ fxseom^ ixk 
■ dc^Mt%' liarticuM^ TO flie descciption and concerts egu^y.i^Iy 

. t» '^iliifef^'con^t^ hefwoi" ^stans, includiag systems having arcfaiteetures dissincbilar to that 

^ ' SfiH r^emng to HG. l.tfe server con^utffl-s^ 

a database server 105. The Web server 103 manages networic resources and handles all 
•/application operations between die browser-based clients 104 and the server side applications, 
TThe database server 105 includes a database management system (DBMS), a collection of 
prograni& tibiM enables the storing, modification and extractioa of informatioii fi'om databases 
114. The Weto SCTV^ 103 ^u^tat^ communication and data exchange between tiie client 104 
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and database 114. 

One or more data acquisition devices 130 can be used to generate:raw^D^ata about an 
object and to input the raw 3D data to the server 102. Some exai^Ies of suitable data 
acquisition devices 130 include laser scanners for acquiring 3D surface data and MRI scanners, 
CAT scanners or Laser Confocal Microscopes for acquiring 3D volume data. It will be 
understood, however, that the data acquisition devices 130 can include any device that generates 
digitized 3D data from an object. 

A program Icemal 116 executes on the server computer S5rstem 102 and implements the 
logical operations of the systan, as described below. The kemal 116 includes a geometric 
modeling module 118, which can derive 3D and 2D niodeled data from the raw 3D data. A 
feature exfraction, analysis and indexing module 120 caii opetafe to pix>vide xatalogingj 
descriptipn and interactive access to the modeled, data, as well as to stored textual and 
.descriptive data about the object. A data compression module 121 can coihpress eeit^n data for 
enhanced storage and transmissioii. 

A region editor program 132 also executes on the server compufer system 102.^ The 
region editor, program 132 ftmctions. to break the mesh down into, ''bbaningful" connected 
subsets of vertices called "regions." : 

The database 1 14 can comprise a nunaber of different datsCbases for storing variotis' data 
element These data elements include: 3D acquired data 122, such as that generated by the data 
acqdsition device 130; modeled data 124, which the geometric modeling module 1 18 generates 
from the 3D acquired data; derived data 125, which the extraction, analysis- and indexing module 
derives from the;modeled data 124, and textual and numeaiq data 126. The database il4 also 
stores metadata 128 for describing the data and how the data is formatted. T|te. sdiema 
..discussed herein further describe, the metadata 128, 

The 3D acquired data 122 can be 3D surface data, such as tbaf generated by opfically 
scanning the surfe|ce of an object, gr it cm be 3I>; yobme ^ta fb^ indlidt^ io&noa^dn' about , 
tfie. iHtCTior of an object, .such as that .Qbtained..fioifi MRI acaas^ 

li^CTOscppe volume cUita, Hie modeled data- 124 is .data tiiat Im b&m stotctered ot indeed 
from &e acquired data 122 u^iug suitable models, sui^- as tdangular m^e^, ^satSice 
representation modds or voIuiae.;rqpre5eatatiOT modeis. llie derived data 125 includes data 
derived from the modeled 63iSst,,si3ch a^- x^jeet, features «t makoiiaficstl descriptive dat^ 
^edsat^ fiffltt fite inbdeled dsta. Zhe .tesat and utmiaEic data 126 can incMe, for example, area 
airv;^e di^tiibntiQn and the like. • . •.. ^ ■ ; ^ 

The data eieruents are orgjunzed, cataloged and deSsait^^ «icik!r#^ 
febilitate efficient qasEypfdMMs^^ 
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the schemas are written in XML (extensible Markup Language). XML is a standard format for 

structured document/data interchange on the Web. Like HTML, an XML"do«uinent holds text 

annotated by tags. Unlike HTML, however, XML allows an unlimited set of tags, each 

indicating not how something should look, but what it means. This characteristic h«lps facilitate 

information sharing. FIG. 2 depicts an example of an XML schema according to the present 

inviention. The schema can include application specific categories and can be orgainzed to 

differentiate general information (e.g., prevalence, time period) and 2D information (e.g., 

height length) from 3D information. The 3D information data structures at the top-level house 

original binary data files, modeled files, and derived features. One level down (lie hierarchy, 

cognitive pattern primitives (CPPs) store the basic shape data that form the building blocks out 

of which meaningful features are constructed. The schema allows shape grammar and shape 

algebra to be used to provide a natural lai^uage syntax to write a shape as a construct of tokens 

(CPPs) and combination rules. To link this spatial^ modeled, and constructed data to descriptive 

data in existing databases, the schema provides a master identification number for use as the key 

to link data elements for a given artifect This model is consistent with Dublin Core cataloging 

structures and enierguig data standards within the discipline. This design permits a query engine 

to link tlie data elements with additional data in other databases by m^^ping appropriate data 

fields between the databases using a master JD#, 

Client Computer System 

Still referring to FIG. 1, the client 104, which also has a conventional architecture, 

includes; a central processing imit (CPU), which may be implemented with a conventional 

microprocessor; means for temporary storage of information, which may be implemented with 

random access memory (RAM)- and means for pennauent storage of information, which may be 

implemeatod with read only memory (ROM); and means for mass storage of infoimation, which 

may be in^lemented by hsrd drive or any other suitable means of mass storage known in the art; 

<me ;,<3(r ifloie k^ut devaees £)r iipittin^ ii^noaliott into ffie cofliputer, wMcli may be 

h^&pasats6 usuig .a "isoiiymtioQ^ k^oaixi aod mo^^ an ouf|)Uf devicd ibr (fisplay 

. jgjtaphical oi^t to tiie user, siich as a coiaventioiial monifbr. Hie cHent 104 includes 

conventionai hacdw?ee aa4 software, &r cdiMur^ali^ the ntetwoifc 106. In a presently 

|irefeEred ranbpdim^^ fivs^^i^^ a Webbiowser sroftware i^)pli<;ation l l0, which executes oa 

&e clien£104 to access Jafoimat^ avaaable on the Memet. Any conventioaal brw IlOis 

cojffitejE^yl^ed for use accor^ng to the pasesart-Mvaatiorii^i^e^^^^ Navigator or Microsoft 

Memet Explorer, provided fbai the ibroTj/iseT^^ and HIML 3.0 or bett^. The 

Systran and method of the present disclosuarc advantageously utilizes the full fimctionality of 

suchbrovyser pr^>g:^ns, as are known topersons skSled in the operation and use thereof 

' • ■ ■ . 11 • ■ 
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Graphical User Interface 

The client 104 utilizes a GUI by which a user can view, store, mAdi^-iofonnatioii in and 
extract information from the database 114. In the presently preferred embodiment, the GUI 
comprises a collection of HTML pages transmitted to the cheat 104 by the Web server 103 and 
displayed by the Web browser 110. The GUI includes interactive surface and volume 
visualization capabilities and quantification tools to extract curvature, volume, scale, and linear 
■ dimensions for each data set and allows the user to generate a sketch-based of the database 114, 

Region Editor 

Referring to FIG. 3, an exemplary cUent screen display 208 is provided firom which 
certain aspects of the region editor 132 may be described- CUent screen display 200 is 
representative of region editor screen layouts that may be viewed by users of clients 104. A 
display window 180 displays modeled data as an image. 182 and numeric data. Tha user can edit 
the modeled data using a toolbar 184. 

Visual Query 

The GUI also includes a visual query interfece,. . which, allqws us^s tp .ii^u^ aiialyze, . 
refine and limit searches of. the database .114. The visual query ioter&ce combines a dceteh- 
based search capability with traditional text and metric data- Kefeiring to FIGs. -4-7, an 
exemplary client screen display 200 is provided from which catain aspects of the visual query 
interface may be described. Client screen di^Iay 200 is represaitalive of visual query screen 
layouts that may be viewed by users of clients 104. The client screen display 200 can display a 
query input screen 202 including, a shape input window 204 and a text ii^ut window 206- The 
shape input window 204 allows a user to .mput a query image .206, which can be in the form of 
vector graphics or a 2D and 3D image. The query image 208 can be created as a fieelbim 
profile sketch or it can be selected from a palette of stored representative object shapes. Using a 
toolbar 210, the user can mampulate the query image 208 m. real time to refine searchesv i The 
text i^ut -window 206 includes i^ut fields 212 for iiaputting textual and numeric search 
parameters. These support parallel query of descriptive and derived data within the databases 
1 14. The user can manipulate and resubmit the query image in real time to refine searches. 
After a user has entered the desu'ed criteria into the querj^ input screen 202j a "submif' button 
213 can be selected to submit the query to database 114. 

Referring to FIG; 5, when /fli© user submits a query; Uie client 104 sends the search 
aiteda to the Web server 103 as binary and textual data (step 250> The system 102 then 
. seaKjhes the databases 114 to match text and mmieric, and cdculated feature data, with the 
seatehisaieaasu- A CoinmoH Gjdev/ay IntaJ&ce (CGQ pmgmm accepts the d^ from the Web 

:l2v ; 
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searver 1U3 and parametenzes the search critena (step 252). The kemal 116 then exti-acts 
features, using appropriate feature extraction algorithms as described bedffsWi-.-fiBmJhe sketch and 
executes a search of flie database 1 14 (step 254). hi a presentlj^-preferred embodiment, Java 
. Server Pages (JSP) are used to implement the dynamic database queries, but it will be 
understood that other suitable technologies such as Active Server Pages (ASP) also can be used. 
The database server 105 queries the database 114 (step 256) and returns the search results to the 
kemal 116. Tlie kemal 116 then compiles and tabulates the search result data (step 258). hi a 
presently-preferred embodiment and method, XML is used represent the search results. Query 
results from the database 114 are stored in XML format, and are visuaUzed via a pre-designed 
Extensible Style Language (XSL) file. The XSL file provides the data needed to view detailed 
information for any of the returned objects without nfeed for further interaction with the server 
102 and databases 114. This includes tllvmibnail images and descriptive data for a predefined 
number of objects in the search results. The Web server 103 sends the processed search results 
to the client 104, where they are cached (step 262) and displayed by client 104 via the Web 
browser 110 (step 264), which extracts fhumibnail and summaiy data Sxim the XSL file (step 
4 14) aiild displays flie tiiambhaii and summaiy dafa ini the initiai search result screen. 

FIG, 6 illustrates an exemplary search result screen display of mitial search results, after 
the client 104 receives the search results, it displays those search results in the result window 
204, The displayed results can include thumbnail images 214 of objects having features 
matching the search criteria as well as and descriptive data 216 for those objects. When a user 
selects one of the returned thumbnail images 214, which prompts the client browser 110 to 
extract a manipulable 3D model of the selected object and detailed descriptive data. The client 
browser 110 then displays the 3D model 218 of the selected object along witli additional 
descriptive data 220 for that object. The user can modify tiie query sketch by retamiug to flie 
initial search result screen or new search screen. To view more details associated with the 
■ selected fhumboail image^ the user can select a detail Enk 221 which will cause a separate 
display page t» be caHed (step 266) and displayed as a full results display 268. FIG, 7 illustrates 
- aa lektaple of a S&resk d^Iay of the fill resiilts" of a search. The fiill r^ults display 268 
iaeludes additionaf 2iy, 3D and descriptive dafa for the selected thutnbiml, including features 
such as pit^ etHfVe data 222, curvatuwiJlot ;data 224 and full textual and numepc data 
226. 

V : : System OperaH^n ^ : 

Referring to FIG. 8, additional aspects of operation of thb system: 100 will now be . 
' des<3ribed in fuither detail. Generally, operation of tiie system includes a, data acquisition 
process 300, modeling of the acquired data 302, segmeartaiion and ifeatiiie enaction 304, feature - 

■ ■ 13- . . .:-Vy 
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anal5'sis and index generation for storage of data 306, and online display and query of the stored 

data 308. 

Data Acquisition _ 

In the data acquisition process 300, 3D raw data is input to the system 100. This raw 

data can be 3D surfece data, such as that generated by optically scanning flie surface of an 

object, or it can be 3D volume data which in<dudes information about the interior of an object, 

such as that obtained from MRI scans, CAT scans or Laser Confocal Microscope volume data. 

In the example of 3D surface data, "vdien an object is scamied i;ising a laser scaxmra:, several tens 

. or hundreds of thousand point coordinates are generated, each representiiig a location on',.the 

object This coUection of points is referred to as a 'point cloud" It has no structure, and is 

simply a file containing point coordinate data as x, y, z values. Point "cloud information also can 

be generated using computer programs, either interactively or automatically from mathematical 

functions, as is well-known in the art. In either case, to work with this point cloud, it lias to be 

structured or modeled. 

Gedmetric Modeling 

The geometric modeling niodule 118 operates on the acquired 3D data to perform 
■ geometric modeling (step 302) of fliat daifa. FIGs. 8, 10 and 11 illustrates tiie data flow and 
structure of the geometric modeKng module 1 18 . Geometric modeling can include mbdeTing of 
3D sutfece data into surface representation models 3 10 as well as modeling of 3D volume data 
into volume representation models 3 12. (FIG. 10) 

Polygonal meshes 

One advantageous method to model 3D data is to triangulate the point cloud to generate 
a mesh of triangles (Step 314) 314 having Uic digitized points as vertices. 3-D triangular meshes 
are a weU'-known surface modeling primitive used extensively to represent real world and 
sgnithetic surfaces in computer graphics. Each triangle 'knovsrs' about its neighbors, v^iiich is the 
structure that allows fast processing of the geometry rqjresented by the triangulafion. FIG. 9- 
iHusMtes an example of a triangulated, hi^y decimated point cloud 6a!t& set 'represaiting a 
surface of an object 316, i.e. a ventral surfece of a stone flake. It is jn^ortant to sfress that the 
same set of vertices or data points could Iiave sev^al triangidations. Tbs axn^Sams ffaetsfore, is 
on the vertices themselves, rather than the 'surfece' as represented by de triangles. 

SwfaceModds 

Because triangdation is a modeling primitive, the ge«Muetric modeling module 118 uses 
additional modeling steps to siBq)iify &e;dat& and extract tenons or object features that are 
msaningfiil fer atppfications or to nseis; For examjde, to stody surfece contour shapes the 
surfecesmHstbe.modeitea. If an object is inliefi^^' ^oolii, ft is represent its 
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surfece by a smooth data fonnat One advantageoiis surface representation model is via a 

parametric surface such an a NURBS (Non-Uniform Rational B-SpKS^-data set 316. A 

NURBS data set coxisists of control points (besides degree and parameterization) arranged in a. 

grid of points, roughly represoiting a given shape. Fitting points of polygonal meshes with least 

squares approximation generates such surface models. NURBS data sets provide a compact data 

representation said m^e it easy to calculate curvatures of objects to extract their essential shape 

•characteristics. Using such r^resentations enables one to rebuild models, arialyze properties 

such as curvatures and msike quantitative measurements, as well "repair" incomplete models. A 

NURBS surfepe can be represented as: * 

_ £ Y^w,.jli.jNi.k(u)Nj.,iv) 

PM^ '-'J"^ (I) 

f-O >-0 . 

where dij, i = 0, 1, . ..,m; j = 0,1,. . .n are control points, mj are weights, Nij/u) aadNjj (v) are B- 
Spliiie basis feuctions. When weights equal 1.0, tiiis reduces to a non-unifonn B-Spline surface. 
2D Geometric modds 

2b models can be used to simplify pioblems presented by 3D analj^is. For exaix^le, 2D 
curves m9.y be used to describe an object such as a profile curve used to describe an 
archeological vessel. To obtain a profile curve, a cutting plane is projected through the 3D mesh 
representiiig flie object, and intersection points with the mesh are extracted and ordered in a 
chain code. NUKB curves can be generated by fitting points of the chain code with least squares 



One of the important characteristics of a curve is the curvature. The curvature can be 
very useful for analysis and classification of object shape. In 3D space tti© ctorvature of curves is 
unsigned. Hbwever,^ for planar curves in 3D space, positive curvatures can be converted into 
signed curvatures. Curvature has useful mformation such as convexity, smoothness, and 
inflection points of the curve, and if this is ioformation-needed for -ana^ 
can be used to j^prosmate profile crarves . , ;. 

m-^„ (2> 

■ ."■.■:/=o,.- 

where di, i ,= 0, 1, .v,,n are control points, ii'/ ^.^wei^ts, and Ni^iiu) are B-Spline basis 
fimcti<>iis. - . In one exempli ^plication, , because measures of ourvatcHre such as- convexity, 
smooflmess, and inflection points of the curve were needed for vrasel .aiMfysis, cubic B-SpKne 
curves have^l>e^ The accuracy of ' 
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these derived curves far exceeds tlie manually sketched or mechanically derived curves currenfly 
used to describe these vessels. The result is a profile curve from which apppopHate- features such 
as comer points, end points, and symmetty between both halves of the_profile curve can be 
readily computed and extracted. 

FIG. 7 illustrates a vessel image 218, its profile curve 222 and the signed curvature plot 
224 of the curve. 

Feature Extraction and Segmentation 
Referring to FIGs, 8 and 11, after modeling of the acquired data (step 302), the feq.^e 
extraction, analysis and index generation module 120 performs feature extraction and 
segmentation (stqj 304) so tiiat appropriate feature data can he extracted for study. Althou^ 
features (such as comer points on a pot, joint surfaces on a bone, and flake negatives on lithics) 
are a diverse set, they can be extracted in accordance with the invention using a common set of 
algorithms. As shown in FIG. 11,. extraction of features from such objects and the use of these 
features to search for particular sh^es involve multiple levels of analysis and extraction 
algorithms. As a common denominator, however, the features are bmlt fiom cognitive pattern 
primitives (CPP), which are essentially the geometrically meaningful feature building blocks; 
Examples of CPPs include curves, for exanqjle a part of a .profile curve, and surfeces such as 
segments of the sur&ce of an object Building features from these piimitives requires a 
construction method that allows flie CPPs to be arranged and connected in various ways and also 
for these connections and arrangements to be describod mathematically, Featares can be 
defined, then, by a s^ of CPPs and a set of operators to connect and arrange these CPPs. 

Extracting features fi-om modeled data can take on various forms depending on the 
desired infonruition. According to one exemplaiy feature extraction method, features are, 
separated into four categories: point, curve^ suifaGe data and volume data. The extraction of 
point features, curve features and surface data features wiU now be, jiescribed \\dfh reference to 
an exemplary application, i.e. the study of G^amic vessels. 

' Extra<M7tgjmm features: co^ 
Using ch^<k>des:aiid]S^^ 
from modeled data, infese ittofile curves can have a nund)ar of points of interest to a user, such 
./s& (6bm^ f^^ Co^t'ip(^xiles cmhs ddsoibed as ah abto^t chai^ m 016 dtieiifation of m 
;i,cbjectwall, or a disfjnet ^le in. the J<Mr^ of object pafts, siK^t\as,neGfc aodbo^ji^of a vessel 
Pot file automatie extraction of snch i»ffits onit'pffiflte corve^^^tf^ is a fl^ed t6 ^Ui^pjafely 
, describe .this.^eatoe m3tfaeni^GaUy/- HeBce^ c(»^ df 'maxiiiimli dOrVatore.- 

dian^ bn a profile ciirve ^ aiti estraisted by coi^jaciag c^B^^^are valne^ fiir ail'pdiiats on tlie 

U 
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profile curve. 

The following algorithms have been developed to extract point fejasfiesr — 

Algorithm for extracting end point features^ 
Input: a profile curve represented by B-Spline curve and chain code respectively. 
Output: end point features 

1. end point 1 □ startpointof chain code; end point 2 □ end point of chain code; 

2. c«iter point □ center point of chain code; 

3 . find the base section around center 

4. //base section is flat or concave then 

total sai point number 0 4; end point 3 □ left terminate of base seetiori; 
end point 4 □ right terminate of base section; 
else {base is convex} 

total end point nmnber □ 3; end point 3 □ center; 

5. calculate feature infortaation for each end points, include space coordinates, 
parameter value, position on the chain code, and so on; 

Algorithm for extracting cortier point features 
Input: a profile curve represented by B-Spiline curve and chain code respectively. 
l9«$?M/; comer poiatfeatures 

1 . calculate curvature value for each points on chain code; 

2. find points with local maximum (minimum) curvature value as candidates for comer 
points; 

3. all candidates 

(f angle at the candidate point < a predefined value 
; tbe caiididate pQiat IS a comer point, 
- 4. cdcuktei^tdieihfbnnatioa breach comer point, kclude space coordinates^ 
'.. . - - ■ parameter value, posifion on the chain code, and so on. 

. Inflection feature^ and point ^of verti(^ tangenf^ features can be determined because 
inflection pou^t and ppp^t.of yertical taagedey can be found by analyzang curvature value and 
tangent lines. ■ ^ ..... , . 

When Qompulirig the angle between points (jq, yO, (xb, yo) aad^ yr) in flie algotiftan for 
: exttac^aig comer features; the aagfe: is smsitive to saa^leerro^ due 
to sampKng, instead of taking (xi, 4814 (xii-.jv) as jpoiiitg.of ^ cmve, th& "ccoidiiates of these 
;poiiits are calculated by avera^g tiie coorditiates^ of a group of m^^Saors to perform a less 
noise-prone resampling. ^ 
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Consider the midpoint (xo, yo) of n contiguous points in a chain code of a curve, where n 

is an odd number, and let p = n/2 + 1 be the point (xo, yo). Thus, the initialpoiat of the angle 

(xi, yi) is calculated from the n/2 + 1 previous point as - _ 




Extracting cuTve features: profile cur\'es. 
Profile curves include various other sources of valuable information. For example, one 
could be interested in the symmetry between left and right side of the profile curve. This 
information can be extracted into a signed curvature plot, and provides a valuable tool for the. 
evaluation of object symmetry. Further an ideal model can be generated from the profile curve 
and used to evaluate how a real object deviates from this raodel. 

Extracting swfa(xfeat%ires: facets oil surf ac^^ 
Another level of feature extraction, deals with surfaces. By extracting surface 
information, one can retrieve 3D inforriiation rather than reducing a 3D phenomenon to 2D 
infonnation. To extract surface information, the complete object has to be segmented into 
discipline specific and geometrically meaningftif sUifeces. Examples of such surfaces are the 
different facets on a lithic artifact such as the ventral surface, partial dorsjfl' surfeces, and heel. 
These surfeees fonn flie building blocks in any effort to refit fiiese' arfifacts into a complex 3D 
jigsa,w. puzzle. EIG.; 13 shows two pieces of an oigect, i.e. a itidte and ffefci? of a. stone artifact,- 
'vSoSii^-^efxs& & t^^ of which pieioes^ hias tii^Bea' lsegpsM^ td 'extmct itis sxixface 

features fer comparison. 

Segpeatafioitin flnsHeoBtext-flais ref^^^td |)E6h&MW or regions 

of -interest &Qm sur&ces. According to one . ail^eet of Ibvention, '&e computerized 
segmeatation of triangular mesh surfices \^&, ?ibis3lufe'6UEvatette estimates along wi^ the basic 
iiieas of a Watershed algorithm to identify re^ns of iiifCTe^ and similarity on the surfaces of 
objeGts. . . Aj^plications of tins method have been found particiilarly useful in recognizing facets 
on a lithic attifeots, sudi as the hthic artifect shown in HG. 13. These regions in turn form the 
basis fot-.seatching formatches between stone aitifacts. 

Based on the idea of a watershed as used in geogra^Ky, m which a watershed forms the 
dividii^ Ene betwe^ drainage basinsy the eo^nteiizea watched segmentation scheme defines 
mininia to. whi<4i wat» would flow.^m fite^pea^ surronnd^ fihaf nmimta, :Ja contrast to a 
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geography application, where this principle uses elevation in relation to neighboring areas as the 

defining characteristic for subdivision, this application uses absolute ^BtKature. allowing the - 

recognition of peaks as well as valleys as the separating lines between watersheds. It was found 

that absolute curvature aided by smoothing equations yielded superior results to Gaussian 

curvature. 

A first step in the se^entation process is the computation and storage of curvature at 
each vertex or data point in the original point cloud. (See FIG. 14) The curvature calculation 
for each point is based on a patch of nine or more points, around a particular vertex and the 
vertex itself. Next, absolute curvature minima are selected and form the "sink-holes" 
individual regions. Subsequently, vertices are assigned to a specific minimum by determining 
the way the unaginary water would travel (dovm the steepest drop in absolute curvature in tliis 
case). The toitial regions are typically too numerous to be usefiil, therefore, we increase the 
"watershed depth" to merge adjacent similar regions (See FIG. 15). Making the action of 
determining the watershed dqjth user defined we ensure flexibility to meet researcher needs. 
While significant depths that allow for the recognition of flake scars are typically found at 
similar watershed depths, these depths tend to be different for other application such as the 
recognition of joint suifaces or parts of a pot such as neck,' body & base (FIG. 15). The 
segmentation algorithm is described in more detail below. 

Metltods for Extracting Curve Features 
Improved Cwvature Estimation for Watershed Segvienicftion of S-Dimejisicmal Meshes 
One advantageous method for ^xtractijjg curve features from 3D meshes accordijig' to the 
invention uses an inq)roved curvature estimation scheme for Watershed segmentation of flie 3p 
.meshes. The melhod improves iipem flie scsheme .proposed in -A. Maogan and R. Whitaker. 
^aititioniag SB Sut^% Medt^ Usu^, W^ateaffi^ Traosacdoiis oo 

Visualization and Computer Graphics. Vol^, No. 4, Oct-Dec 1999 (Mangan and Wiitafcer) for 
segmentiiig a 3D |)QlygoiiaI loesh r^^a%senta£iQa of m arbitEmly shifted real vr&M object 

The domaia of the.piobleei^^sHeKi a in£sj& (d^ 
(verdc^ Vt e B^;. 0 ^ i < a) and a sef of plaimr convex polygoiB (faces) made up of these 
vertices. 3D . meshes iare currently a popular surface modeling primifive, used extensively to 
represent real world and synthetic surfeces in compute gtapMcs. In oMer to gMerate ineshes,^ 
real world objects are often either digitizedj ie. point samples m:e collected frbin ihe surface of 
the object and then ccHmectivity generated or sampled as a volume fixDin which aa isosurface is 
extracted. Alternately, Ak points (and fheir coimectivity) are generated using computer 
programs, either interactively (e.g. using a surfece desiga tod) or automafically (e;g, fi-om a 
mLathemalical function). A mesh is a distal or. ^iiu'^re^zzed' structure r^seseating some : - 
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continuous suiface. 

Segmentation means breaking down an existing structure into jrreaningfiil, connected 
sub-components. In the context of 3D meshes, the sub-components of the structure being 
broken down are sets of vertices called regions which share some "coramonaUty." A value l-, 
could be associated with each vertex v,- in the data set which somehow encapsulates the 
diaracteristics of the locality of the vatex. The dej&nition of segmentation is one in which 
regions consist of connected vertices which have the same X (within a tolerance). Curvature is 
selected as the mathematical basis for region separation, i.e. the scalar value X. Previously, 
curvature estimation &om 3D meshes has been dealt with by extracting curvature from a lodally 
fitted quadratic polyaomial approximaut or. by various discrete curvature approximation 
schemes. Vertices having, tlie same curvature value would be grouped into regions, separated by 
vertices with high curvature (which serve as region boundaries). 

The segmentation scheme described herein is derived from flie watershed algoritiun for 
3D meshes described in Mangan aid ."Whitaker. The watershed algorithm for image 
segmentation is well known in the 2D iinage segmentation field. Mangan .and Whitaker 
generalize the watershed method to arbitrary meshes by using either the discrete Gaussiacn 
curvature or the norm of covaiiauce of adjacent triangle notnials at eadi me^ vertex as the 
/lejgAr field, analogous to the pixel intensity on an image grid which drives the 2D watershed 
segmentation algorithm. Altfaou^ the npyel algorithm described herein is derived fromMangan 
and Whitaker, it differs in flie curvature estimation scheme detiails and is far more stable for real 
world noisy data. 

Surface Curvature 

The quality of results from the watershed segmentation algorithm depends substantially 
•on the accuracy and stability, of the-«stimated curvature values. Curvature at the mesh vertices 
should faithftilly reflect the local properties of the underlying surface. Additionally, \X\b useful 
to have a . curvature calculation techm<pe that can a) vesfea^ eurvaiur© aecuEately oa rdesh 
boundaries where there usually aren't enough v^ces distribtifed evettly around the vertex in 
question, and b) is rdativdy unaffected by noise in the data. 

To, provide an .understanding of the in5)roved curvature .estimation scheme, varigus terms 
used in liie tteory of surface fiBTvaturewiU bo briefly desaibe^ 

.Ilifi^r<niiefric srarfeces consid€ied are of the form 

jc=;<u>;u=(tt,v)e[a,b]cR^ (1) 

^ere » and v are parametes vi^udi take real values and vary freely in the donmia [a, b}. 
The fimctions x(i^ v) = <x(% v), y^, v}^ z{u, v)) are single valued and continuous, and are 
assumed to possess continuous partial dmvatives^ 
20 
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The first fimdametital form, denoted by / is given by 

I^x'x^Edv^^lFdudv^Gd^ (2) 
where - 

^ = x^ = x„.x„, F=x„jr„ G = ;cv = JC,-Cv- (3) 
The second fmudamental form, denoted hyUhs given by 
n^Ldu^'+lMdudv+Ndv-, (4) 
where 

X = Nx^,M=Nx„„7V=Nx„. . (5> / 

and N is.the surface normal at point X. 

The nortnal curvature of the surface at point x in the direction of tangent t is given by 

Kn=Kn(x:t)= = — (6) 

Since the normal curvature is based on directioii, it attains maximum and mimmmn 
. valuesj, called the principal curvatures. Hie principal curvatures, Ki and YLit. can bfe combined to 
give us GiiMSJKiB curvature, givai by 

Gaussian curvature is invariant under arbifeaiy transfonnations of the (m, v) parameters of 
a surface as long as the Jacobian of the (m, v) transformation are always non-zero. Gaussian and 
mean curvature are invariant to arbitrary rotations and translations of a staface becauise of the 
invariance of the E, F, G, L, M, N functions to rotations and translations. 
Discrete Curvature Schemes 
_ . . For the purpose of comparison, we first consider the two discrete curvature estimation 
schemes used in Mangau and Whitaker, i.e. the discrete Gaussian curvature scheme and the 
. flOHn of file, coysffiance of the sa^we noimals of taaagles adjacent to the vertex at which 
curvaticeishei^ealeulaled, ' ... 

: pjs,crete Gaussiaq ...... 

; ; . The platelet of a vertffl£V< is defined as the Sdt ofaU tnangles iii-tilie t^esh sharing Vf as 
a vertCK- The discr^e Gax^ 

(8) 

where is tte angle subtended at v,- by eadi tiian^^ey in P/, and^l/ is the area of that feiangi©. 

A problem fliat presents itself inaaediately is, how to fflid theplatelet of a y^ex on the 
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boundary of a mesh that is not closed (e.g., a mesh representing a sphere would be closed, i.e., 

have no boundary, while one representing a bounded plane would not be ^dosed). . The platelet of 

a boundary vertex would be incomplete, giving rise to an inaccurate computed value for 

curvature. One option is to simply ignore boundary vertices while computing curvature (setting 

their curvature to some fixed value). However, this oft^ gives rise to unpleasant results during 

segmentation because such vertices tend to establish their own regions distinct from the rest of 

the mesh. 

A solution to the problan would be to extend the mesh beyond the boundary in some 
"meaningfiil" way, allowing a curvature computation using the extended platelet. Suda, an 
extension is not trivial. ■ It would have to somehow follow the shape of the mesh (and 
additionally, it is desirable that ttie extension be symmetric) in ordar.for the boundary vertices to 
blend correctly into the existing mesh, 

Gaussian curvature is m isometric invariant property of a surface. An isometric 
invariant is a surface property that depends only on the E, F, G fimctions (and possibly their 
derivatives). Isometric invariants are also known as intrinsic surfiice properties. An intrinsic 
proiperty does iiof care how the surfeoe is embedded in higher dimensional space. On the other 
hand ofiier curvature metrics such as mean and absolute curvature are an extrinsic property 
which does depend on the onbedding in 3D space. Inliinsic piopaties do not change sign when 
the direction of ttie normal vector of the surface is reversed. This is because the first 
fundamental form does not depend on the surface normal, while te second fundamental form 
does. Hence, Gaussian curvature maintains its sign when the difectiou of the normal vector is 
flip|)^ whereas mean curvature (described below) flips its sign. 

Isometric invariance is an undesirable property as far se^entation is concerned. For 
. example, in FIG. 16, we would ideally want the .left and rigjit.halves of the sujtface segmented 
info separate regions at the ridge; The Gaosffli^ curvature plot sho#s that flie curvature is 
. constant thrQu^out. The segmenfatfOn, being based on ctirVature (K^ifributioi^ would be unable 
.^- tq.dKtiijguish the halves; Th^ mean -ciaviEiture plot, on the other hah^ sho^ro the ridge being 
: /distinBtly mai|;ea by Mgfrcittvatar&== iy an ad^tioiffi^ disadvaiitage, Gsiissiii curvature being 
the product o£&e|£ri^^ 

Norm Method ;•• '•• '' ' ■ . / 

Hie second method used in the original p^o: computes the. noim of fii0 ;CWariflnce of 
the surfece normals of the trian^es adjacmt fo the vested at wj^eh cwrvai^ is b^iig calculated. 
■Rris inefliic^ we sJ^/^p^ to as ifee Noim lai^od, fe awre lesffieait to noise in ^e data. 
It is bet^ #an QSmg Qansdan bti£ does not perfbim as wdl metrics used in ihe improved 

- ■ ■ ■ 22- - ^. : ■ ■ ■ . ■ . ., 
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curvature described berein. This also has no direct geometric significaace. 

Lnproved Curvature Schemes 

The principal curvatures can be combined in other usefiil and geometrically meaningful 

ways such as mean, RMS and absolute curvatures. 

Mean cwvature is given by 

~ 2 2 £-G-F' 
The mean curvature, being the average of the principal curvatures, is less saisitive to 
noise in numerical computation than the principal curvatures. . 

Root mean square curvature (RMS) is. a good measure of surface flatness and is given by 



2 

and can easily be computed as 



k^=4aH^-2K ill) 
A value of Kthk = 0 a* a point mdicates a perfectly flat smface in te 

poiiit. 

.<4£>J(3/Mfe curvature is given by 

K,^=hil + N2| (12) 

Mean and RMS curvatures have a closed form as given by equations (9) and (11) and 
these do not require actual con^tafion of principal curvatures, ko and k\ are expensive to 
compute and . hence tiie cost of computing Absolute curvature is considerably higher than 
computing mean and RMS counterparts; 

Curvature Estimation 

Curvature estimation from a triangulated mesh can be a non-trivial pffoblem depending 
on flie accuracy required and the method used. Our method uses a technique •wfiiere the surface 
is approxunated Ioeall)( by a biquadratic polynomial whose second order mixed partial 
derivatives are calculated. These in turn aUow computation of the E, F,- G^ L, M, N fimctions. 
An added advantage of fitting a surface locally is that flie sensitivity of flie surfece to noise can 
be "acyjusfed" by using smoottung equations. 

The pareunetric sur&ceusedis liiensor product Beziersiaface and is^giveaby: 

:c(«-v) = 2; i;h,X(«)5;(v);«,v€[0,l] (13) 

bjj are called Bezier control povits: in their natural ordering they a*e the vertices of the 
B62ier control net. -m:.... 
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We perform a standard least squares fit to a set of vertices in order to compute the bij. 

The following were used to improve the local fit of the surface (FIG. 1S)^~--— - 

1. Use double star of the vertex as input data points. _ 

2. Add smoothing equations to the least squares solution 

. Selecting more vertices around the vertex for which the curvature is bemg computed 
ensures there are no gaps and the resulting surface does not behave wildly. The second option 
invokes constraints that a good control net would satisfy. One such constraint is to nunimize the 
twist of the control net, which, for a surface is its mixed partial ^/dudv. The twist surface of 
b"^" is a Bezier surfece of degree «-l) and its vector coefficients have the form mnA^^bij, 
■ Geometrically /BnA'"'by measures the deviation of each sub-quadrilaferal of the Bezier net fixjm 
a parallelogram. Minitnizing the twist has the effect of making each sub-quadrilaterial close to a 
parallelogram. Hie Constraint can be stated as follows: 

A'-'b,_^. =0; 0<i<m-l, 0<7<«-l. (14) 
From Equation 14, the condition can be restated as 

Finally, fitting a surface of hi^er degree than 2 gives the surfece . greater jfreedom to 
match tile underlying data points^ However,- suifeces of M^r degree come at a higher 
computational cost Moreover, fit offered by. any sur&ce of greater degree than a quadratic is 
usually a case of diminishing returns. A second-degree surface suffices for evaluation of second 
order partial derivatives. 

The first tiiree solutions presented above miake the Imear. system qyer-determmed. This 
system can be solved using least-squares methods. HowevK, we need additional information 
i.6. parameter values associated to each vatex. This is solved as follows. 

The data points (vertices) can be parameterized by projecting them onto, a plane and 
scaling them to the [0, 1] range. If the data points (in 2D) . are distributed.,^ shown in FIG. 17A, 
vvith the mia-max box as shown, scaling them to the [0, 1] range results in areas in the domain 
wifb. no data points in &em. This is undesirable as explained previously. To get a better 
parameterization we obtain as small an enclosing box as possible. This can bedone by selecting 
a local coordinate frame such (hat the min-max box is miaimizsed. This problem can be solved 
by finding the axes of the smallest ellipse that completely encloses Ifie points. Ihe axes of such 
. m: ellipse, called the norm ellipse, would then result in tiie smallest min-max box for the point 
set(s©sJlG,17B). 

The Watershed Algoriton 
; The watershed algorithm witt now briefly Be 
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: in Mangan and Whitaker and the literature cited therein. Once the curvature ki at each vertex v,- 

in the mesh has been computed and stored, the segmentation can begia:^ -KrheFe represents a ■ 

generic curvature metric and can be Gaussian, absolute, etc. During segnentatiori, a label will 

be assigned to each vertex v2 in the mesh indicating to which region the vertex belongs. 

The 3D watershed segmentation scheme aims to segment the mesh into regions of 
relatively consistent curvature having high curvature boundary vertices. It is instructive to 
imagine the curvature distribution "ftmetion" (or flxe height function) as a surface on the given 
surface. The watershed algoriflim operates on this (curvature) surface. 

Convex areas (elevations) on the object have a high magnitude of curvature with a 
positive sign, wMle concave areas (d^ressions) have a high magnitude of curvature wifli a 
negative sign (except for curvatures which are positive by definition). High positive curvature 
would appear as a ridge on the curvature distribution flinction,. while a high negative, curvature 
would manifest itself as a valley. Since the original watershed segmentation algorithm segments 
the sur&ce into regions separated by areas of high curvature, the valleys would not act as region 
separators as would be expected (see FIG. 19A). We only considered the magnitude of 
curvature (FIG. 19B) to solve this problem. 

Results and Conclusions 

A flexible segmentation program was implemented, with user control provided for 
changing parameters and thresholds in order to arrive at the desired segmentation of a surfece. 
Here we show the results of our segmentation process. 

Among continuous methods, mean cvuvature was more resilient to noise in numerical 
computation siace it averages the principal cutvattires. Being an extrinsic property, mean 
curvature produced results closer to the expected results which were based on user perception of 
shape, RMS curvatare had its advantages whea deafing with specialized segmentation. 
HoweveF, absolute curvature resulted in sdj^entations &at goieratly outperformed all others. 
Like mean, the absolute curvature is tiie smaaa&m of (Jibe absolute values of) iq and j<2 giviog it 
greater noise resiS^ce. Moreover^ it is positive by definition unlike mean and Gaussian 
cuFvatnresi. 

GMMparison of Qgvatnres and their Est imatian Techmg nes 
• •Continuous curvatare-esttetion methods resulted & more accurate cuivature estimates 
than discp^ mettiDdSi This was because: 

♦ Surfaces could be fitted to largCT vertex ndghborhoodg, increaang the closeness 
of the fitted surfax^ to the actaal undfiriying surfece. 

• Smoottung equations tended to offeet the effect of missTTig data points and high 
fceqaencyixMseia^e data, set 
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« Boundary vertices did not need special consideration as in the discrete case since 

the surface could be fitted to more vertices on one side of the hoandary. 

Segmentation using different curvature metrics. Use of Mean. and_RMS results in almost 
the same quality segmentation as use of Absolute at a cheaper computation cost. 

Segmentation of human trq)ezium (tiiumb.joint) was performed using three methods for 
curvahire approximation. Discrete Gaussian and Norm methods resulted, m too many regions 
inconsistent with the geometry of the surfece. 

Segmentation of a lithic using three methods for curvature approximation. Discrete 
Gaussian and Norm mefliods produce geometrically inconsistent region where as absolute 
curvature results in correct segmentation. 



Worst 


Gaussian (K) 




. Mean (H) 






Best 


Absolute {Kabs) 



Taible I : Comparison of curvature measures for segmentatiGn. ' 
From the foregoing, it can be seen that an advantageous method for segmenting a 3D 
mesh has been presented using a 3D mesh watershed segmentation scheme. The segmentation is 
based on curvature, and therefore several curvature estimation techniques have been 
implemented and analyzed. The robustness of segmentation has been improved by increasing 
the accuracy of the curvature estimates by use of mean, RMS and ahsolute cutvaturesii Tlie 
method successfully segmented several real world data sets, both digitized arid syBthe^e, .and 
often ?ioisy, primarily because of the eflfort that went into developing better ctirvature ^timatioB . 

. iie m^od jpffli be implemeiited in ajapgram is. semi-automatic since fhe-rfeSult is 
(Iq^dent CKQ IJie mbr^t^ threshold uspr iopiit :Op&,pf,&et reasoa&Aidiy.a stu^ fiif^^M'<!o^ 
iwt on.all.dat^ sets.:is:€ia£ ;sinSc^ caa^ in.varioM siz^i and; 

curvature is no^t sp^g ittvaiiant .UptfoEmly. aalin^ homMis^hfm of the input 

surfeces to a fixed size aUeyiates flie pro ft was aoted fljat sorfec© 

sets of .die same tj^e 0^g., ppte) i»spoi^wi|ed the satHettoeshold 

;^ues. The other reason for this is. tlia4.i .pdfectiy :¥a^^ disagree wiflj die 

•%spected" user perceived , segtaentation. Hiis was especially, tmf in this research- since the 
cementation goal for data by the requireiiierts of 

experts from theses disciplines. V . • . :• - - ■ - - . . : 
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A Hybrid Approach to Featwe Segmentation of Triangle Meshes 

Another advantageous method for providing feature segmentation" ef- triangle meshes 
according to the invention will now be described. Previous approaches for segmentation of 
polygonal meshes use a vertex-based method. A significant drawback of the vertex-based 
meliiod is that no hard boimdaries are created for the features or regions. Each A'ertex of an 
object has its own region infoimalion. Therefore, triangles on boundaries have multi-region 
information. The tiiree vertices of a triangle can be part of three different regions, whereas the 
triangle itself would be a "gray" area, i.e. it woidd not belong to say one region, FIG. 20 shows 
an example where a triangle T/ will be such a gray area. FIG. 2GA shows an exarhplb of 
seginentation using the watershed method with no hard.boimdMes. In FIG. 20A the boimdary 
triangles are shown in the white region. This means the regions will not have hard boundaries or 
edges. This is a known artifact of vertex-based Watershed segmentation and is acknowledged in 
Mangan and Whitaker. To provide a solution for this "no hard boundary" problem, the vertex- 
based method of feature segmentation creates new triangles by adding mid-points to the edges 
that have different region labels. FIG. 2GB shows the segmentation exiample of FIG. 20A- using 
the watershed method with such a boxmdary solution. 

FIG. 21 illustrates the creation of triangles on boundaries of an original mesh according 
tQ the hybrid method. FIG. 2IA shows the creation of trisoigles on boundaries of an original 
mesh for a triangle shared by two regions. FIG. 2 IB shows the creation of triangles on 
bovindaries of an original mesh for a triangle shared by three regions. Referring to FIG. 2 1 , each 
new vertex contains multiple labels. This is an extension of the original algorithm. For the 
selection of the label(s) of a vertex which has multiple labels, the eomihon label of the vertices 
of the triangle is selected. In FIG. 21 A, there are two possible diagonal splits to make triangles. 
We selieet the (fiagonal whit* satisfies max^tn\^i£4/au 7WOTfi<i4bx)f where a^aad ^^ are interior 
an^es Ifonned by the diagonals^ Le. select the diagonai that rssiilts in the best aspect ratio. 
.i|i^S0e vds^es jDqjreseidiBg^echsueai (CAI^ bt>]^^a»& jlecpt^tty s^e^ &i Some isess, 
^fli just eoot^ y^ees to define GsaSk^e&; Tbifiis i^soally a rgsulf'df tfie optnaial tiiangtilation 
fioni « CAD progtam or d^dmafion'piDdcess. M- this case, Wiit<a:s&ed me^od may not 
-se^eitt the objjsct prop«ly or may even lose Hiq)dttant rej^oiis on objects. Bi HG, 22, die 
.r^^ tffeated iKbdTimdariesbiecaase tiiere are not enough vertices m this 

re^ensi- The boundary-solution mentioned above will iiot soiye this problem because llie method 
n6f creafe n^ re^ons .fiom the houndaiy regioifs' aii^ lie hoimdarj^ i^ons vwll be lost, 
MoreoYCT, some regions of this object are not segmented properly. This problem is caused by 
ihe v^ces on featine edges wifh ^tmlar curvatntes and those v^Ces may fee treated as the 
J same re^dn. We solve Oris problem wi& our piK^ 
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This method uses an edge-based method for defining bonndari^iS:r. .ArEeature Edge is 
defined as . an edge shared by two planes whose nonnal vectors make an angle greater than a 
certain threshold. The edges obtained are integrated into curves, and these'curves are classified , 
as jump boundaries, folds (roof edges) and ridge lines. Jump boundaries and folds are used to 
segment the mesh into several regions. The boundary lines are also treated as Feature Edges. 

The main disadvantage of the Feature Edge-based method is that this results in many 
. discoimected Feature Edges and thereby incomplete Feature Loops. FIG. 23 shows this problem. 
Feature Edges are shown in the shaded spots, • ' 

Hybrid Approach 

The hybrid method will now be described^ This method creafeff regions with complete 
Feature Loops. 

FIG. 24 shows the steps of the Hybrid me&od. As explained Before, optimally 
triangulated meshes pose problems for the Watershed segmentation mefiiod; To overcome this, 
all of the feature edges in the mesh are identified using a threshold (angle). First, a F.eature 
Vertices are defined as VCTti,ces that make up a Feature Edge. The vevexse not-n^icessarily true; 
■ if both vertices of an edge are Featuie Vertices, it does not automatically qaitijfy the edge as a 
Feature Edge. 

FIG. 25 illustrates the creation of triangles on feature edges of an original mesh. 
Step 1 - All Feature Vertices are first defined based on a threshold angle (dihedral angle) 
as explained above. 

, Step 2 - Next, new vertices are added. Vertices are added to edges of all triangles that 
have all three vertices as Feature Vertic^. See FIG. 25. The new vertex is added at the mid 
point of each edge. We then coimeet, than as shovsfn; to. create four new triangleis. -If the triangle 
ha? dree Feature Edges, then the center point of the triangle is added and six new triangles are 
pleated 3s shpiOTi toPIG* 25B. / Addiiea of yertiees reqtiiress fij^ ©f fiie topdogy ^ fflusMtei 
in -^G. ?5Cv.^Tnie.tiiaD^? skovmis ai a^ijhboF of ttte triage to wMdi we ktdedi&e v^ces. 
.Ti^s x5j^ i^. to a banging pioMem. To fix -fiui^ eoimiset ftk'i&^ ydiek'-)l^:1}i6 . 
■YWt^ .on.the.ppppsite edge to .create two. new triaBgles: As a; result of iflie above, we hkVe two 
Jiiads^oriiew ve£ti<^^^^^ fhat donof 

^P^ied as FVIowX The reason fertabdib^^ ::■ 
, . Sie . ; 

, Next, wQ ,^ply Watershed segm^eettfltti. ttf o«r modffied mesh. We yssQ ^tbsdute 
Guwai^re for our exargil^^ O^er curvatuie estiaatioa SKSiods may be ijsed Feature Vfertices 
FV Kgh are ass^^-.fiug label of laajdniani <^imtoe. Since Ihey lie on a Feature Edge, 
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assigning them high curvatui-e ensures that a Feature Edge will be preserved as a hard edge. The 

rest of the vertices in the mesh follow the same procedure as described^ -the Watershed 

algorithm for computing cmrvatures at the vertices. The Feature Vertices contain their own 

region labels as well as labels of the neighboring vertices. The addition of vertices has an impact 

on the Descent and Region Merge operations of the Watershed process [9]. This is done to solve 

the \no hard, boundary" problem, wbich has been described and discussed before in this paper. 

Step 4 - Removing Vertices Added in Step 2 . 

To restore the mesh, to its oijginal form we must remove the vertices added to the mesh 
in st^ 2 above. This process restores ttie topology of the mesh also. '.' 
Step 5 - Collating Triangles into Regions 

The goal of this step is to assign triangles and not vertices' to different regions. Tliis is 
adiieved as foDows: 

1 . Case: All vejrtic^ have the same label 

This is siiiiplest of the cases. The triangle is assigned the region label of its vertices, 

2. Case: Qiie vertex has a single label . 

This is the case when one vertex has a unique label but flie otii^r two vertices, have 
multiple labels„.Tlj.e .triangie is assigned the re^on of the vertex wi& single labet 

3. Case; Multiple labels but only one common, label 

The three vertices of the triangle have multiple labels eachj however, there is only one 
label that is common. The triangle is assigned the region label fliat is common to the vertices in 
the triangle. 

4. Case: All edges are Feature Edges 

The triangle qualifies as a region by itself and gets assigned a unique region id^tifier. 

5. Case: Multiple labels and multiple Ciwnmon labels 

. It is possible that the each vertex of a triangle has multiple labels and there is more than 
oae common labeL To explain this we will use the exaa^e in HO. 26: Trian^e Tl has veitiees; 
with region labels Rl and R2. .0ne of the v^ces has Ifie laBel R3. First, we. check if a 
nei^sboiii^ triangle shares u Featisre edge with HI ; la Qd& case^ we find T2 shares common 
. Feature edge with Tl . We them coropare ^e. oonmion vertex labels df Tl (Rl R2> with coinnwn 
label(s>^)f T2 (Rl>. The label that does not Jbelong to the set of common vertex iabel(s) of the 
nagywiing triangle is assigned to the tjffgeted tri^ 

If rfie targi^ed triangle has more flian one feature e^ 
for each nei^boring.triangte . that shares a feature edge -2 . Each such niMgK)Or <X)a&ibutes one 
or more, labels. Thea.Tl is assigned the labetthat is ©ommoa betwe^ the contributed laifaels.^ 
The targeted triangle may have no featare edga It means that the tSaiigle is in the saine region 
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as the regions of the neighboring triangles. SOj the region label of any neighboring triangle is 

selected as the label of the targeted triangle. 2 Note that if all three edges-ace- feature edges, then 

it is dealt under Case 4. _ 
Results and Conclusions 

The hybrid method relies on the Watershed method and additionally uses advantages of 
the dihedral angle method. The hybrid method does segmentation of smootli objects as well as 
mechanical objects. It also solves the "no hard boundaiy" problem. FIG. 27 shove's two 
mechanical parts. FIG. 27A is a gear and FIG. 27B is a housing. FIG. 28 A shov/s a turbine. 
Notice that each blade has been segmented into a different region. FIG. 28B is a wheel (rim^pnd 
tire). It is segmented into the following regions i.e. the tire, rim, bolt holes and a hole for 
mounting the wheel. These are correetly segmented. All exanq)les are from pubKe domain data 
sets. 

HG. 29 shows the bone data, (from FIG. 23) segmented using the hybrid inethdd. We 
have pointed out the shoitepmings of the most popular segmentafion approaches fer triangle 
meshes and devised a hybrid method tliat overcomes these shortcomHigs. Futwe wo]& Tvould 
involve automatic selection of threshold values both for watershed and dihedral angleis. 
_ Data Qompression Algorithms 

As previously described, according to one aspect of the invention, tlte data con^ression 
.module 121 corapresses modeled dak for enhanced storage and transmisision. Two 
advantageous compression lE^thods are (1) triangle mesh compression ustBg B-spHne ' curves 
apd (2) compression usuig a; subdivision surface scheme. Eaeh of these compression methods 
will now be discussed. 

Trmngle Mesh Compression Using B-Spline Curves 
As more sophisticated objects are modded -or digitized,- their -triangle mesh . 
representation becomes more complex and hence the size of such data increases considerably. 
Due to their large size, these files are difficult to transfer and tafcp a long time even for viewing, 
-Triangle Mesh Compression, techniques exploit , thq geometric properties of the mesh 
representation to encode them into smaller files.. A novellossless compression scheme using B- 
Splines is described. 

A triangle mesh represent^on consists of infcnmation about geometry and connectivity, 
alsQ known as /(^(j/ogy. Geometry defines the location, of v^tices in a (Euclidean) coordinate 
system. They are represaated as triplets (x^ y, z). CoimectiyLtj' defines the sets of ^joints tiiatare 
tsounec^sed to form triangles or feces of the mesh. Itianglea sffie-^vm bjr.-ftEee.indas vsflues, . 
which identify the three vertices bounding the trian^e. 

(jeometiy information is H5ua% cosgaessed by qoaQti?ii^ the vertices. Based on the v 
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number of bits chosen to represent the; floating point numbers, tlie location of the vertex is 

converted to integers on a grid. Note that ev^en though there is loss in the" gSGmstriG- accuracy due - 

to quantization, it is still referred to as Lossless compression. Lossless refers to the connectivity 

i.e. the triangle mesh is reproduced after decompression topologically exactly the same as the • 

input mesh. A suitable prediction scheme is used and the displacement between the predicted 

and actual point is stored as an. error vector. When the prediction is accurate, the error vectors 

are very small and hence can be stored using fewer bits. Connectivity compression is achieved 

by designing a traversal scheme, vAibh encodes the unique path tafeen during inesh traversal, 

xtshig predefined rules and tables: ,< 

The present novel conq)ressiou sdieme vises' B-SpUnes and is Idosely based on the ideas 

•presented in the Edgehreatoer algorithm described in. J. Rossigpac^Edgebreaket-TraossiCtiojis on . 

Visualization and Computer Graphics, 1999; J. Hossignac and A. Szymezak, triangle meslies 

compressed with Edgebreaker, Computational Geometry, Theory and Applications, 14(1/3X 

1 19-135, Novanber 1999; and J. Rossignac, A. Saifonova, and A. Szym, Edgebreaker on a 

Corner-Table, Conference, Gatnoa, Italy. May 2G01 (Edgebreaker). A brief descriptioiilof the 

Edgebreaker algorithmfoHows. 

A "seed triangle" is randomly picked fiom- the miesh. Starting from it, the mesh is 

traversed one triangle at a time, based on rules of traversaL The edge through which the seed 

triangle is entered is called a "gate." The path taken is written as a combination of codes. When 

a triangle is entered, it is marked as visited; all of its vertices are also marked as visited As We 

move from one triangle to another, the ferther vertex of the entered triangle is predicted based on 

the three vertices of the first triangle, using parallelogram prediction (FIG. 30). The difference 

between the predicted vertex and the actual vatex is then stored as an error vector. 

, , ; The error veGtors. are stored only for vertices wMeh have.not been visited yet. Hentie^ 

after the traversal, fhe^ number of codes wotdd bb; equal to tie number of triangles and ttteie is 

exactly oae exrorveGtOfeozcsspondiog to each v^teK; ' ' 

, The topology is raided as a s^es of codes coirespooding to the mesh traversal, and the 

gemnetry is encoded as error vectors. Usnig tirist ia&mBs&ott, ihe original mesh is neconstruoted 

during decompression. 

DiHingomr tests wenoticed lhat topology encoding m up for roughly 10% of the 

eompr^ed file (coa^jieesed nsing-Edgebreaker). The chaUenge th^i is to reduce the 

coaq)r!essed file size Telafed to ihe geometry: Haioe tisae-is a need for a better prediction 

algor^iiDi than the parallelo^Eaat iagaiod. We ejqperaaettted vritfe B-Spliaes. "Hie duid benefit of 

using our method are tiiat even though we use the mBsih traversal method of Ed^reidcer during . 

compresaon we do not stare miy'codes^ iind, &e enor^ptediction is improved The eonnectivity 
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is the automatic result of the algorithm. First we describe the compressidn. 

CQinpression :^r-^ 

The B-Spline compression scheme uses the same traversal rules as described in 

Edgebreaker. However in our method, when we move from one triangle to the neighboring one, 

the mid point of the common edge is recorded. Once all the triangles are visited, many 

sequences of triangles, called triangle strips, are obtained. Each such strip has a corresponding 

sequmce of edge mid points, to which a B-Spline curve is fit using least squares approximation. 

The order in which the E-Spline curves are geinerated is recorded (FIG. 31). The parameter 

value, corresponding to each midpoiut on the curve, and the error vector, which is the difiference 

betwe«i tile actual position of the midpoint aid the one computed using flie parameter value, are 

also recorded. 

The parameter value is stored only if the vertex opposite to the gate 4as not been visited 
yet. If it is marked as visited due to the traversal of an adjoining; fmii^gle «arliOT. in the . 
compression process, the Judex of the edge that is shared with that atSpiniBig triangle is stored 
. for referehdng if diirihg decompression. This avoids storing the already wisi^ed vdiices multiple 
■times. 

Decompression 

After the cdmpressi^n process, we obtain an ordered set of B-:%Eiiifi control poiiits, the 
parameter values corresponding to- flie midpoints of edges and the eaicEtr-'VBCtor associated with 
each mic^oinl. The two v«tices foiming tiie gate edge are also given. 

During decompression, the B-Spliae curves are evaluated at theiir!E£Etrf parametisr values 
in the order ifi- which tiiey were created. Each pokit thts evaluated on literiJitBRe is liieB displaced 
by the corresponding error vector to generate a mid point (FIG. 32). TihisiDiai poin^ along with 
one of the two poinfe<jbtaiaed dur&.g the previous sf ep (iiOring tiib &s3t^^ 4hey would be tte 
gate vertices) is used^ deterttune tlve tyri vertex fori^ 

. Points that are revisited store a refereitce to the edge commion ftnsfeesa the current 
. triangle and the adjoinmg one-visited previously. Tfiis e<!ge is 'used ttwlleKsfllB fc? a<^ual poin]^ 
whicfa is-theii used to construct &e cuir«it triangle. 

Preliminary results obtained by coa^ressioii using the B-^jifiBB sdv^aie are very 
. mcoiHaging. lie control polygon points and ffie enor veclois are stojE^ -as tmhs of floating 
point numbers, each stored iJsSrig 8 lafs. The huaibe^ sStoiBgEaraai^ter vdues is 

dynamically computed based on number of p oints on the B--Spline cin?«eBianafe typca^ 5 bits. 
Tte B-Spline . curves used for our experimerrfs are of degree f&«e sand 1^ 
parameterizafion(FIGs. 33anid34), , 

We note the foUomig two idfflferteac^ 
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about both geometry and topology is encoded into the B -Spline curves. Hence, unlike 

Edgebreaker, our method does not store codes for topology. Second, the ^axoi^iuectors associated 

with mid-points using B-Splines are about ten times smaller than the error vectors obtained in 

Edgebreaker, which uses parallelogram prediction. On an average, the compressed file size is 

10% smaller than that of Edgebreaker. 

Compression Using A Subdivision Surface Scheme. 

Subdivision surfaces are finding their way into many Computer Aided Design and 
Animation packages. Popular choices include Loop, Catmull-Clark, Doo-Sabin etc. Subdivi^sion 
sxirfaces have many design advantages over traditional useof NDRBs. NURB surfaces always 
are problematic when multiple patches meet. 

Reverse engineering (RE) is associated with the idea of scanning physical objects and 
representing the resulting dense, cloud of points wiflt mathematical surfaces. In RE the goal is to 
convert the dense point of scanned points into a patchwork of NURB surfaces with most effort 
going into automatiiig the process. With &e emergence of subdivision - surfaces as popular 
modeling tools, it follows that a similar process be devised for this class of surfeces. Our 
p^er looks at devefeing one such method for Loop surfeces. 

Given a triangular mesh, we would like to.obtain a control inesh for a LOop 

subdivision surfa^Pphic^ approjomates the giyeii .mesh. This process benefits subdivision 
suifiices in aninia^Q and manqjulation that need speed over accuracy with an . ahilit^^ to 
manipulate the 99^1^1 mosh. and to regenerate the smpofli sgr^e quickly. ISks subdivision 
wavelets, in a mul®esolution analysis, owr mefiiod. can perform leyel-of-detailg, (LOD) with 
arbitrary topologic^ meshes useful in applications requiring a fiisit transfer, less storage, and a 
rmdedng aid interaction. The reswlt^ coarse control njesh is approxi^nately 6.25% of &e 
m^i^ rosi^ therefore this method can, also be used as a lossy con^jr^saon scheme. 

Giv«i dense miDrgani2M d poii^„,sQcli a?: a poiiit. cloud firom:a range scaimer, a 
trifflQgle mesh can :|e constructed by various methods known in. the git However, triangular 
meshes obfeined ar$ piecewise linear surfeces. For editing, modelings etc. dense triar^lis; meshes 
are not an optimal solution. . Dense triangle messes with tot of detail are expensive to represent, 
store,: transmit and manipulate. A tensor produ^ct Ntj^ Rational B-Splines) 

and B-splines are die most popular smooth sur&ce ipprcs^itation^. Considerable work has been 
done on fitting B^Kne surfeces to ftBceerdhnQBsioml points. This process is often called the 
,5Lc?Varsefiigineeriag process Tviwaieina^talj^iresentation of the jiiysical obje<*is tasated. 

A B -spline surface is a paranietiic siBfecev^ Previously, 
some have proposed methods for pac^meterLdng a triangular mesh and- im(a:]^iiized points 
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respectively for a single surface patch. A single B-spIine patch can only, moael suriaces -mtST 
simple topological types such as deformed planar regions, cylinders, ^d^tOTi, Therefore, it is 
impossible to use a single non-degenerate B-spline to model general closed surfaces or surfaces 
with handles. Multiple B-spline patches are needed for arbitrary topological surfaces; however, 
there are some geometric continuity conditions that must be met for adjacent patches. Therefore, . 
using NURBS/B-splines is not the most desirable approach since it requires high-dimensional 
constrained optimization. A subdivision surface scheme such as Loops, does not suffer from 
these problems, has compact storage and simple representation (as triangle base mesh) and can 
be evaluated on the fly to any resolution with ease. It does not require computation of a dopiain 
sur&ee. 

Similar in natufe to subdivision wavelets, we would like -to obtain a control mesh 
approximating original surface with small magnitude of details (as a scalar ftmction) needed to 
beist reconstniGt an approxiination of the original mesh (FIG. 35). The benefit of using a scalar- 
valued function is that its representation is more compact than its traditional counterpart, a 
vector-valued geometry representation. Oiie of our contributions is to obtain a control mesh with 
y&y small magnitude of displaced values (details) in order to have a very high compression ratio 
,while preserving surface details given a semi-regular (subdivision connectivity) of the original 
mesh. 

Subdivision surfaces 

A subdivision surfece is defiiied by a refinement of an initial control mesh. In the limit of 
the refinran^t process, a smoolb. surfiiice is obtained. Subdivision schemes previously hkve been 
introduced based on quaifcilateral meshes. These schemes generahzed bi-quadratic and bi-cubic 
tensor product B-splines. The triangular based subdivision scheme was introduced by Loop 
■which was a generaliaation of quartic triangular B-splines, 

Remeshing 

Remeshing means to convert a general triangle mesh and into a semi-regular mesh. It is 
called a serai-regular mesh because most vatices are regular (valaice six) exc^t at some 
vertices (vertices txm a control mesh not having valence six). 

Semi-regular mesh can be used for multi-^iesolution analysis. Through the remeshing 
process a parameterization can be constmcted and used in dfiier contexts sucfi as texture 
m^ing or NURBS patches. ^ 

The. overview of our jnefliod isvas followisz an otiginal mesh is sdbt^rlified by a quadiic 
error mehie (QEM) based on fi^oii^ 

enxff metrics. Proceedings of . 24th. itramal aaaSBcaasst aa. Cioo^puter gr^hics and interactive 
techniicpes August 1997 to get a simplified control mesh. We dedmate tiie inesh to 625% of the 

■ ■■■;.34 
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number of original triangles. Based on Suzuki, H., Takeuchi, S., and Kunai, T. Subdivision Surface 

Fitting to a Range of Points. Computer Graphics and Applications, Prdcefedings. -Seventh Pacific 

Conference, 1999, the control mesh is then adjusted such that after some levels of subdivision the 

control mesh vertices he close to the original surface. The adjusted control mesh is then 

subdivided using the Loop's scheme Loop, C. Smooth subdivision surfaces based on triangles. 

Master's thesis. University of Utah, Department of mathematics, 1987 for two levels obtaining 

subdivided mesh with its number of triangles close to that of the original mesh. Based on Lee, A., 

Moreton, H., Hoppe, H. Displaced Subdivision Surfaces. Proceedings of the annual conference on 

Computer , graphics (SIGGRAPH 2000), 2000, the subdivided mesh vertices are then displac^ to 

the original surface. Our remeshing process is shown in FIG. 36. Most or all vertices, however, 

require displacements during this process even if the control mesh has been adjusted previously. 

With that in mind, we need to reverse the process of Loop subdi vision scheme to obtain a better 

control mesh with smaller and fewer displaeement values (details). A flow chart of our process 

is shown in FIG. 37- 

In the Loop subdivisionj .each, subdividing level produces two types of vertices: vertex 
points and edge points as shown in FIG. 38. After obtaimsg the subdivided and displaced inesb, 
the mes3i is semi-regular, (subdivision connectivity) and all. vertices Ee on the oiigincal sitrface. 
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visavcrtcxpoMt- 
« is an edge point . 
I is subdivision level . 
Afis& Loop subdivision matrix 



The first reverse step applies to themedi to obtain a coarser 1*' reverse subdivision mesh. 
We don't choose a subdivision matrix in FIG. 39 A because both vertex points and edge points 
are used in solving linear system of equation. The subdivision matrix is not a square matrix, and 
we need to do least-square approximation via the nonxial equations to obtain the contrpl^verticeg. 
If fliat has heea done, d&placed values would have beai required for all vertices (bQlb vertex 
points and edge points), and by doing that, it defeats our. goal of having fewer nijniber of required 
detail values (displaced values and error vjdues) and smaH magiHtjrfQS needed f^^ 
cornpressioiL Therefore, oidy the stihdivisiiaB mabm toe vertex points as shovra in FICj. 39B is 
tised instead, Here, the matrix is a square in^rix: said; is inyertibte assiHning.,a.Yprtex. gene^'al. 
position. W& cgkM solve farexai&f cootrnt vertiees; For krge%^i^^^i^ii^^^^ 
iteratiye ^raacli to sglve ffie^ Enear Systran. Naw, we need fe'<Mcd:8rie-^s0ai3ed-vaI||^ and 
save fhent for a reconstrucfion phase. To gettheoii we infernally subdivide flie I^Ubvms© 
subdivi^on meshvoae- levet The vettex points of (Ms subdivided mesh (levd 1) are about 25% 
of total vertices, and they all have zero displaced values. The edge pomts may or may not 
require displaced values. The percentage of edge points requires zero displaced values as shown 
in Table 2. Note that the coBJ^Jutation for displaced values for tedge poii^ at this level is to find 
the distance along each vatex limit normal intersecfing vrith the (o^inal surface. 
M the second revise stqp, we agffln solve for &e conbol Vffliices of f^^ 
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in similar maimer to that of the first reverse step. However, the error values (or what we call the 

displaced values in the first reverse step) for the edge points are compttteardi£fer«itly as shown 

in FIG. 38. The error values are computed as a dot product of edge pointy unit limit-normal and 

vector &om its position (obtained fcom subdividing the solved eontrol medi) to its 

corresponding position (obtained firom the 1'* reverse subdi\asion mesh). The end result is a very 

sr&all coarse mesh plus some details using wiiidi one can approximate the input mesh. The 

coarse mesh can be used as a model much like the NURBS/B-spline control mesh or the model 

and. detail can be used as a lossy compression scheme. 

We reconstruct and compare a number of well-known data sets as shown in FIGi' 40. 

After obtaining control meshes, we losslessly compress them using software by 3D Compression 

Technologies Inc. to further reduce the ou^ut file size. For the encoding of details (error values 

and displaced values) we vary the quantization level fi^om 8 to. 12 bits per de^ value for 

comparison. Magnitude of Loop vertex limit normal is used as a scaling &ctor to flulher reduce 

the already small detail values. We compare the result between 8-bit and 12.-bit detail encoding 

and found that by encoding with lower bits (8 bits), the reconstructed surface is as good as that 

of the hi^er bits (12 bits). In additfdh, 8-bit quanti2ation would produce higher number, zero 

details, \vhich can be further encoded using variable-length scheme. In the variable-length 

encoding, each detail value is either encoded by 1 bit for zero detail value or 9 bits otherwise: 

The ninth bit is to indicate the required detail encoding. As a result, our output sur&ces have 

high fidelity of original surface details as well as high compression ratio. Table 1 below shows 

•percentage of vertices rcquiriiig no displaccdmcnt values at different levels. Comparison ol: the 

quantitative compression results among 8-bit detail encoding, 12-bit detail encoding, and 9-bit 



variable-length detail encodmg COPTZ) are shown in Table 3. 



Data 


Level 




% of no> of vertex 

points require zero 
displacement) 


%«f oo.af edge 
poktstf^tbzero 
displacement 


Total % of no. 
Of points 
requiring zero 
displacement 


Spock* 


Control mesh . 


.1,039 




16389 




4,121 


25.2t 


15.97 


41.18 




. 2 


16,417 


^ . 25.10 


.62.65 ■ 


87.75 




^Confrolmest' 


2.102 




. 3^87 


. i . 


8.398 


25.02 


8.86 


33.88 




2 


33,586 


25.00 


20.47 


45.47 . 


Biuiny+ 


1-ControI'mesfai' 


2,206 




35,280 


- X 


Mis : 


25.02 


13.60 


38.62 : 






35,266 


2500 


24.34 . 


49J4 


Horse* 
48,485 


CoKbxiI mesh 


;:3;032. 




1 


ixm 


25,01 


55.02 


. 80-03 




2 


48.482 


25-00 


72.15 


: 97.15 



Table 2: Zero <SspIacemeat at each level (widi 8-bit quantization). 
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Data 


#V 


Size 


Control mesh 


Details (displaceme'nlj 


X^mpress 




#T 


(KB) 


#T 


3DCP 
(KB) 


after 
3DCP 
(KB) 


(KB) 






(control 
mesh.+ 
OPTZ) 




32,718 




1,039 
2,044 


36.13 


.3.54 


15.87 


22.87 


6.21 : 


98J0% 




33,587 
67'l70 


1,181 


2,102 
4,198 


73.82 


5-85. 


31.18 


46.18 


27.15 


97.21% 


Bunay + 


35.280 
70,556 


>1,24.0 


2,206 
4,408 


77.5 


••■6,1:6::;; 


32.5 


48.5 


.26.76 


9735% 


Horse* 


48,485 
96,966 




3,032 
6,060 


106.5 


I'm"::: 


44.5 


.67.5 


:^5125V. 


99.09% 



Table 3: Quanti^ve compresaiqa results. 

Legend: 

Data with * courtesy of Cyberware; + courtesy of Stanford Univ. Computer Giaphics Lab 

#V = Number of vertices 
#T = Number of triangles 
KB = Kilo-Bytes 

3DCP = Lossless compression by 3DCompress.com sofUrare ^erfonned on control- mesh.es) 
OPTZ = Optinu2Bd eacodmg scheine (vanable-lehgfh encoding) 

According to the novel method disctosed; we can approximate any arbitrary topological 
boundary and hon-bo\indaiy surfeces with high compression and details. We have combined 
subdivision sur&ce and scalar-valued displacement for our sur&ce reconstruction. Our domain 
suriace is obtaiaed in such a way that it is close to the original surface, the magnitude of 
di^laced values t^ds to be very smaH and is &vor^Ie as a conqaression scheme. With the hig^h 
compression and &ithfiiUy d^Wled sva&pe, the imethod is can be used for applications such as 
mesh compression^ animatioi^ suxfiice editing and man^ulation. 

CtiheMcai^Uig Algorism 
According to another aspect of the inv^tion, the systan can perform curve matching to 
(^iopatB the siaiiilOTty or dis^sin^anty of two shapes.. As used herein "sh^e" means a two or 
tfiree-ditnessibnai' curve. One of the in^)ortant characteristics of a curve is its curvature. 
CJurvatuie iS very -u^eM for arialysis and classification of curve shape. We assume A and B to 
be free fohn curves sudi that, eidier cinvatures can be cosapvdei, or these can be approximated 
by a fimction or a paramestric curyfe of dei^ee tWo of hi^^. juSD ^sce the curvature of curves 
; is unsigned. However^ for planar curves in 3D space- we can r€5piesent curvature as a signed 

■ quatatify. ■ 

■ ■ ' A curve niafchmg ii^dd acc^^ con^ares two curves to obtain a 

■ mfeasuri "wiii(i4ecides the similarity on^^t^^ The ; 
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curvature plot of the curve with smaller chord length is slid along the curvature plot of the 
longer one. The difference in areas of the curvature plots is obtained.- -and >the minimum 
difference value is attributed to the amount of match betweeii the curves. _ 

The method of curve matching has been implemented on the profile curves of pots, 
which profile curves are constructed by fitting a least squares curve to the points obtained by 
intersecting a 3D vessel with a cutting plane. 



The curvature of a curve measures its failure to be in a straight line. The faster the first 
derivative X'(t) turns along the curve X(t), the larger is its curvature. Curvature is defined aS flie 
rate of change of the unit tangent vector X'(t). For a straight line, the curvature is zero as its 
second derivative vanishes. For a circle, llie curvature is constant. Let X(t) = (x(t), y(tX z(t)) be 
a curve, then the curvature at t, k(t), is defined as 



where X = dX/dt and X = d?X/dt^ :and dsBOtes tiEie. cross prodticL A three- 
dimaisional. cmve has non-negative curvature. For planar curves, we may assign the signed 
ciffvature^, 



The points at which the curvature changes sign are called inflection points. Therefore, a 
planar curve may have inflection points but, a three dimensional curve has none. 
Curve Matching Algorithm 

Given two curves as input, the curve-matching method according to the invention 
computes a measure that determines the similarity between the curves. For planar curves in 3D 
space, the cmvature plot of the curve Mvajg the smaller dior(iI«ig;& is slid alottg the curvature 
plot of the otbor curve. The dt^reace in anas of the corvaturs plot 'is compoted by an 
integration routine sucfi as the Trapezoidal rule or Simpsons rule. The mitiimiim difference in 
area is noted and this, along witii a weired component of the difference in chordlengths of tJie 
curves, produces a measiue wMch detenoines &e similarity between: cUrv 

Given two ounces A and B to match, the first step is to ev 
curves, lie two chordlengths are eon^aredaid flie minimiEtn of the two is used to select the 
lengfii of file unitintervaL What this means is that, the two craves are discretised into inten^ais 
of the same len^ The selection of the vmt interval length deteimines^&e precision of the 
result The gteala the number of discrete int€iTaIs,tiiemca«acciirate^^^ 



Curvature 
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where, is the chordlength of curve A, Ig is the chordlength of curve B and N is ihs 
number of discreet intervals. 

The parameter values of the two curves are then recomputed with respect to their 
chordlengths. To find the new parameter value, the following procedure is used: 

A graph is plotted with the chordlength as the X-axis and the parameter value as the Y- 
axis. At discrete intervals of the curve, the corresponding parameter value is computed by 
interpolatiag the con-esponding graph pomts. This procedure ensures that the curvature pipt of 
the curve is drawn with respect to the chordlength. 

After obtaining the new parameter values of the curves, the first derivative and second 
derivative are computed at the discretised points. They are used to compute the curvature at 
these points. These are the curvature function /a aMfg of the curves A and B respectively. 

The curvature plots of the two curves are then analyzed to measure the similarity 
betweai the two curves. The minimum difference of the two curvature fimctions determines the 
measure for curve matching. 



This can. be viewed as sliding (this is discreet sliding) the curvature plot of the cun'^e with 
smaller arC length alorig the curvature plot of the other, and the difference in curvature plot areas 
is calculated at all positions. Since, the curves and their curvature plots are in discrete form, 
with ftie same index representation, this is a matter of finding . the area by trapezoidal rule. 
Different methods of computmg this difference in area can be used, such as the trapezoidal rule, 
Simpsons rule, and Riemann sunt The tre^zoidal is the simplest For better results a more 
approf>riaie merhod hue/ be used. The BuninnHn curratnre area differeaee (dreamiu) md its 
positioa (a) are noted fixim the a^Kjye process. 

Let (lA,fA ) and (/j, j(i the leag^ and curvature fimctt&ns of th©,twe ciarves^ and B. 
Then, the total differmee<<p.wfll be, ■ 




(4) 




In the above formula yj is a factor contributing to the difforence in chordlesngfhs of tlie 
curves and ^2 is the factor contributing to their shapes (curvature), A wei^ted combinatiQn of 
sh^e and size gives a good metnq to match fee curves. A percaitage match value (pemc^c/?) 
can be conqjuted using the following formula: 
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It can be observed that the percentage match is always between 0 and 100. The^ 
percentage match can be used as an important factor when querying-:a::Tiatabase for similar 
shapes. The best case occurs when a curve is compared with itself. In that case, the difference 
in the chordlengths is zero and the difference in the curvature plot area is zero. The value of the 
metric is zero. This is a perfect match: A = B. In. other cases, when a curve A is compared with 
another curve £ having a different chordlength and curvature plot, the difiference in the 
chordlengths is positive and the difference in the areas of the curvature plot is also positive. The 
value of the riietric is positive. The metric determines the amount of similarity. 

The different combinations of tiie weights si and J2 in the above equation leads to 
diffwent types of matches, i,e. partial matches, exact matches and shape matches. Variations in 
si and S2 and each of these diffCTent types of matches will now he described. 

The first type cmve matching that will be described is the partial match. Given two 
curves A and there is always a possibility that one curve is a part of die other. But where 
exactly dp they fit each other? Sliding one curve over flie other is quite computationally 
intensive because the orientation of the curves may differ. By using the meflipd of curvature 
plots, one can slide tbs curvature plot of one curve over the curvature plot of the other curve, 
which is similar to sliding one curve over the other. Since the method uses curvature, which is 
transformation invariant, the probloBs caused by orientation of curves are avoided. For a partial 
match of curves, the shape and size of the curves are considered. In a partial match, tlie 
chordlengths of the two curves may or may not be equal. The curvature plot of the smaller 
curve is slid along the curvature plot of the bigger one. The difference in area is computed using 
the tr^ezoidal rule, bi a partial match no weight is attributed to the difference in chordlengths. 
The minimum difference m area and its corresponding position are noted This is the best 
possible partial matdbu A threshold can be taken as a u^er input, wMck detemiines how close 
tibe curves should be. lliie nietric for curve similarity is as follows: 



where d is the measure for curve matehii^ and^ and /b are the curvature fiinctioHS of the two 
curves^ and .B, respectively. If Qievalueof the^easute^f ^0, liieji.<4— 5. If the value of flie 
measure <f </&-as/i<7/^, then c: 5. 

The next ^e ofcurve matching that will be discussed is the exact match. In the case of 
partial matches if was observed that the curvature plots and tibieir andysis ate a hetta: tool to 
ddennine the siinilarity between curves. However, in partial matches, fte diffetetiiGes in &e 
cjiordlengtiis of flie two curves was not .tak«i into cansideiation. . In the exact match type of 
niatclHng both tfae^^b^ and the size of the curves are taken into consideration while 




(7) 
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determining the similarity between them. The curvature plot sliding is used as the main tool 

here to determine similarity, hence, making this method transformation iuyariaat.- 

In the exact match of curves, the chordlengths of the two curves may or may not be 

equal. Therefore, flie curvature plot of the smaller curve is shd along the curvature plot of the 

bigger curve and the minimum difference in area is calculated by the trapezoid rule and the 

corresponding position noted. This position is the best possible partial match. A weighted 

difference of the chordlengths is added to obtain a metric for curve similarity as follows: 

. d = ^sl(l^-I,y + S2minf''y,-fBt (8) 

where, d is tlie measure for curve matching, is the weight attributed to the difference in 
chordllengths, S2 is the weight attributed to the difference in the areas of the curvatare plots, Ia is 
the length of curve .4,7^ is the length of the curve J3, and are the curvature fimctioiis of the 
two curves.;^ and .5 respectively. The value of the ineasme d determines fitesiiBilfflifr|f between , 
the curves. 

The third type of curve matching is the shape match. Matching two curves by thdr 
sliape must be size independent For example, one may want to = compare the similarity between 
a smaller circle and a bigger one. In this case, tiie shape matehing metriG should:-ah.ow a perfect 
shape match. Thus, to perform a shape match, the input curves are prerscaled. to fiie same size. 
The scaling is done witt respect to the chordleagh of the corves. The process of scaling is as 
follows: 

Given two curves A. and 5 as inpute,^ the first step is to compute Ijieir QhLordlengttis. Let 
U be flie choKUeDgth of the first curve and let be tire chordleE^ of Hie second curve. Let s be 
the scaling factor. Mafhemsdcally, j is defined as 




(9) 



The curves are initially discretisei and Ihe oarvafure values are computed' at the 
jifisGretised points. The curve with smaller ckordlength is now scaled by the~scaling toctor s, 
which essentially means that the chordlength -ralues of the smaller curve are scaled by s and the 
. curvature values are scaled fay l/y. After the process of scaling is completed, the difference in' 
areas of the curvature plots is con^puted by one of the integration-routines previously discussed 
This idififeence in area is metric detmnining the amiladty between the two curves. The 
metric for curve similarity is as fixlkkws: 

d^mm^[fA-fBf (10) 

if the value of the measure if - OjBaea curve A and ctffve £ are of jamilar shape, litis. 
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two similar curves with different sizes be compared for shape using tliis method. 

The foregoing method obtains a measure for matching two cUTve§^~Thaniethod provides 

an estimation of the shape similarity by the property of curvature of each of the cnrves. The 

difference in areas of the curvature plots of the t^vo curves is estimated and this measure is used 

to determine the similarity. If the value of the measure is small then the shapes are similar and if 

it is large, they are dissimilar. According to the invention, this curve matching method 

described above can be used to query data sets for matching curves. 

Example Uses 

The foUowiag examples serve to illustrate certain presently preferred embodiments' and 
aspects of the invention. Ttese examples are not to be construed as limiting the scope of the 
invention. 

Three of the examples describe projects that involve archaeological and biological 
material, namely ceramic vessels, lithic artifeefs, and bones. These projects are: {I) "3D 
Morphology of Ceramic Vessels" which has be^ used for acchiving and searching information 
about the structure of ceramic vessels; (2) 'Xithic Refitting" which has been used to (partially) 
aiAoniate the cefLttsig process ^uougjh 3D scanning and surface modeling; 2ai3. finally (3) "3D 
Topography of Joint Surges" which has been used to automate segmentation' of osteological 
features and quantify the surface area and cun'^ature of joint surfaces and the congruency 
-between reciprocal joint sur&ces, allowing for the diovelopment of biomechanical models. 
Coinmon to all of these projects are the following aspects of the invention: geometric modeliiig, 
feature recognition, and the develbpioi^t of a database structure -aimed at making 3D models of 
these iffdfects available online for query. 

Example 1 -3D Morphology of Ceramic Vessels 
3D knowledge plays an important lole in archaeology. For exffliq)fe, archaeologists 
study the 3D foim of ..Native Amaican .pottery to draracterae tiae dev^pn^ of cultures. 
Convention^y, vessel classificatioa is doae by an expert and is subjective pr&ne to 
inaccuracies. The Hieasui:emeiits are crude and in some case eyeballing is the m^od-of choice. 

.-. As part of the 3DE, project at Arizona St£^ University in Tempe, Ari2X»na, a system 
~ according to the present invention has been developed as a pilot project for the study of the 3D 
morphology of prehistoric Native . American ceramic vessels. The aim of this pilot project is to 
learn about vessd uniformity, and variation with respeetlo fenction as indieatois of developing 
craft specialization and complex social organization among prehistoric Native American 
cultures. ;To accurately c^Uire the complex curvatures and vessd form and size variation, the 
use of mefeic reffiias and Visual jnq)ection of 2D profiles is inadcsiiBate. Hiae&ie, fer the study 
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of prehistoric pottery traditions, scalable, visual, and quantitative comp^sons of .curvatures 

have been developed according to the invention. . 

Archaeological vessels were scanned and defined as a set of three-dimensional 
triangulated meshes composed of points, edges and triangles. The original data was then 
modeled with parametric surfaces, extracting features to raise the level of abstraction of data and 
organiziag vessel data based on XML schema. A Web-based visual query interface permits 
-.users to sketch or select sample vessel shapes to augment text and metric search criteria to 
retrieve original and modeled data, and interactive 2D and 3D models. 

Key identifying features of the vessels were identified to support seareh and .d^ata 
retrieval, and a catalog of terms was developed to describe these features within the context of 
anthropological description, cataloging standards, and emerging " metaidata classifications. 
Mathematical models were developed to translate these features into measurable descriptions. 
. -Software was developed to extract features from the vessel data and to raise theJevel of 
abstraction of data. An additional result is geperationof vessel measurements fer more accurate 
thjm has been possible using the traditional tools of anthropology. 

Shape information is obtained fiom scanned tfare©- dimensional, data of the 
ardiaeological vessels,- usii® 2D and 3D gecmetric models to represeiit scamied Vessels, 
extracting feature information firom geometric models and. storing the feature ioforiiiatibn'in a 
database for Web-based retxievaL A Web-based Visual Query Eaterfece (V QI) is used to archive 
and search vessels. 

Based on 3D scaimed data, 2D measurements (height, rim diameter, minimum diameter, 
maximran diameter) and 3D measurements (area, volmnej vcdume of w^l) were generated, and 
a measure of symmetry was. developed. The pilot project used a control collection consisting of 
two Mexican Jlower pots, two Mexican mold-made pinata pots, and; three .haBd-i)uilt 
Tarahumara jais to verify the validity of the formulas generated. The prehistoric collections 
studied consist. of a total of 87. Native American vessels fi-om the Tonto Basin, Aiia^naj and 
Casas Grandea in northern Mexico. The control study and the two sets of prehistoric vessels are 
used to generate a quantitative analysis of vessel uniformity and symmetry. The' automated 
measurements are important ; improvements* which facihtate the quantitative assessment of 
ceramic unifonnity sad standardisation, inificators of fije devdopmeat of craft specialization, 
differentiation of labor, aij4 tt® developing of ccnnplex ixas of social orgaaizatioii among 
preMstoric culturess. . . . 

Fordiita iacguiation, each vessel was scanned usi^ a laser scanmr to^ 
a given vessel The data acqiuation devipe 13.0 mcluded two seamiKB. ISie fiistis a €yberware 
Model 15 scimner, which is amolHlescaana: having .a-TO<^ .<^afe0Bt.75 xM smd is 
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useful for smaller objects. The second laser scanner is a Cyberware Model 3030RGB/MS 
scanner, which is large stationary scanner allowing the capture of large 43bjex:ts.-£oth of these 
. scanners are marketed by Cyberware, Inc. of Monterey, California. Both scanners are equipped 
with a turntable on which the object to be scanned is. mounted. Resolution and accm^acy of the 
scans are a result of the distance between the scanning head and the object The large scanner is 
therefore less accurate than the small scanner. Specifics on accuracy arid resolution of the laser 
scanners are avail^le from Cyberware, Inc. While individual scans take only 17 seconds, the 
total average scanning time for each vessel is about two hours, dq)enduig on complexity, color, 
texture, etc. The scanning produces data as a highly dense "point cloud" of informationvthat 
visually describes, but does not physicdly represent the vessel. 

The vessel data was modeled with parametric surfeces by ovefliaying a three-dimensional 
triangulated mesh onto the point cloud data to define the vessel as a set of three-dimensional 
triangulated meshes composed of points and triangles. Post processing included filling of holes 
(due to scanner *'oversi^t"), removing noise, etc. A BaH Pivoting algorithm was impl«neQted 
to convert point cloud data sets into triafigle roeshes; In some cases light decimatioji wais 
performed to reduce the density witihout losing accuracy of the overall structure. The result is 
valid point set and topology data for each scanned vessel. The resuK is a model of the vessel 
that is composed of parametric sur&ces wifli physical, measurable tittributes.^ 

Mostly archaeological vessels are (approximately) sor&ces of revolution, and studying 
contour shape will suffice to gather shs^e infonnation about fSae whole object, Aixording to 
archaeological definiticai, there are four kinds of feature points on profile curves to calculate 
dimensions and proportions of vessels. They are End Points, Points of Vertical tahgehey, 
loflectioii Points and Comer Points found on the vertical profile curve of a vessel. End Points 
^s) are points at the lim (lip) or ^ the base (Le. top and bottom of vessels). Points of Vertical 
Tmgmcy (VTs) are points at the place Vfh&e is the maxrnmni diameter on spheroidal ferin or 
Tnimmnm dtameter on bypeibolic t&m. In^tbh Points (IPs) are pomts of ohaiii(g|€i'fiom 
.. ceneaiVe ta convex, or vice versaL €omer Points- (CPs) are poiitts^ b^ a profile 

curve. FIG. 43 illustrates these firar fcmds of f^ftiEe paints of vessel profile euryea 

Four features are cQumm to all Vtesels, L6. drifi^&j mx, hoiy Base. Otiflce is the 
<ip€&Dg qf fite yesseii or fl® miniTntim ^diatnefer of tt© ojferaag, ^^Ms^mii^'he the same as the 
. iin^ -OF luelow &e ifei. Rim is the fbisbed edge of fee top br ppMag of the vessel. It may or 
may not be fee same as fee orifice. It may have a larger fBame^. Body is fee form of fee 
Vessel heiow fee oijfice and abbve- fee hasp. ':■ Base is iS^ '^ii^jiM oi ^ y^sei, portion upon 
iwfaich it rests, or sits <^ a sus&ce. Tfiist niay 6e c<mye^ ffcat, er concave, or a cbi&binadon 
of^tese. HG. 44,iHnsteates^^fe i©.r ^>»nvex t>asB; (b) flal bj&e wife zero 

. . " AS 
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curvature; (c) concave base; and (d) composite base. 

From the foregoing definition for characteristic points and ccffinnQa. features for all 
vessels, feature representation of vessels was formalized as follows: 

<Point FeaturO :=<End Point FeaturO |< Point of Vertical Tahgeacy Feature> |< 
Inflection Point Featiire> I <Comer Point Feature>; 

<Curve Feature> := <Rim Curve Feature> j< Orifice Curve Feature> |< Base Curve 
"Feature>; 

<Eim Curve FeaturO :=<End Poiiit Featurex End Point Feature>; 

<Orifice Curve Featiire> :=<Comer Point Feature> <Comer Point Feature>; ^ 

<Base Curve FeaturO := <End Point FeaturO |<End Point Featurex End Point 

Feature> <Region Featm'O :=< Neck Region Feature> |< Body Region FeaturO (< Base Region 

FeaturO; 

<N6ck Region FeaturO := <Rim Curve Feature><Orifice Curve lea 
<Body Region FeaturO := <Orifiee Curve FeatureX Base Curve Feature; 
<Base Region FeaturO :==:.<®aseCiirve FeaturO;. ■ ' 

<VoIiinje FeaturO :=< Unrestricted Volume FeaturO t< Restricted Volume FeattorO. 
-XML was used, to represent information of vessels. An XML schema^ as shown in FIG. 
45, was designed to rq)resent geometric information, feature information and measured Value of 
archaeological vessels. Feature, information is extracted fixtm. geometric informa^on and is 
organized accorcHng to the feature formalism in the XML sehana. Also feature inforinatton. is 
used to index vessels stored in the database. The XML schema fixr a sample- vessel witii values 
isshowninAppaiidix 1. 

As discussed above, the result of the 3D laser scan of^tbe v^el and initial 'processing is 
. .a.polygonal mesh comp.Qsed of feceSj, edges aod yotices and4hjeir connectioHs. Surface models 
..are generated from the scatfered poiiits of this polygonal mesh by -least sqttares fitting and/or 
.. rotating profile curves. B-S^line suri^es are used to represent &ese sui^ice models. The B-. . 
Spline models are used for model rebuilding, error analysis, 'elosirig holes and !g^s fiwRMi on the 
archaeological artifects, and measured value getting. 

The mathematical modeling algorithms described berein were used to pass surfaces 
through scanned point eloud data to generate measurable data &om- these relatively small, 
diverse data sets. Surface and volume modeling algoritbniS were ^^plied to modfii md generate 
quantitative, descriptive data about the artlfect The data: and processes developed grew firom 
and aie consistCTt wiflLthfi descr^five vocabujaiy of cecamiestese^h^asin i^^fl^pology. The 
level of accujaqr in documsitation and measurematit of the Jfffifects fe exceeded traditional 
techniques used in fte fidd., Hie faittay-.and deiived data about the Vessel provide a record, that . 
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can be re-analyzed in the future and provide a tool for research without physical access, at 

remote locations, or after repatriation of the vessel. FIG. 41 shows tcvjl -examples of polygon 

meshes with watershed defined areas. FIG 41A shows a complete vessel and FIG. 41B shows a 

partial-vessel, FIG. 42 illustrates feature segmentation of a complex shaped vessel. 

Archaeologists analyze vessels by defining various regions using a profile. The profile 

curve is obtained by passing a verticsd plane perpendicular to the base through the vessel. 

Typically profile curves are sketched free hand, often by tracing and duplicating half of the 

vessel to create a symmetric, but not necessarily accurate, representation. In this research the 

polygonal mesh model is used to goaerate a much more accurate profile curve than has been 

previously possible. The resulting profile curve is processed to remove noiise due to scamiing 

error , at the lim. Vessels are initially cataloged into four broad shape categories - simple, 

composite, inflected, and complex. 

A segmentation schemes based on curvature, as described herein, was used and therefore 

several curvature estimation techniques were used. The robustness of segmentation has been 

improved by increasing the accuracy of the oiirvatiire estimates. Due to the accuracy required 

for. good segm^ts^on, computing curvature is a coinplex, non-trivial tasL In this rese^arch, 

multiple curvature estiniadon schemes were used- and compared. The estimation schemes can be 

broadly classified into two categories — discrete and continuous. Discrete Schemes extract 

curvature directly from the geometry of the mesh by estimating the local deviation of the mesh 

from. a flat surface. Continuous schemes first approximate- the mesh vertices locally with a 

polynomial surface, allowing calculation of various forms of curvature from the resulting 

analytic surface by methods of differential geometry. 

The ori^nal scanned data, and the modeled and calculated data have beeij linked to 

existing jecord sets containing the traditional descriptive data about location, provenance, and 

collection of the vessels. The XML schema uMiig a metadata schema derived from the Council 

fm the Freservatioa of Anthropological R&cotis (€OPAR>: ■ Was used to caOalosg aiEHi b^ganizc? 

2D and 3D vessel ddta:11i6 schema defies ^emetifs for a gjivoa si&&et, andli^ data 

aoross miuitiple databases was developed^^ind in^lemenled to st^port res^rcH ipleries. The 2D 

data &om exis^g dat^ases can. be iiicorppiSted by u^ng the scheoia to transit and link the 

search processes wia databases housing the 3D data. T3ib project provides access to data' from 

fee e:idstmg ceramic vessel datab^^ ^rionai ispalial and volume data acquired for 

.scasaed vessels. 

To^provide efficient access, a scanned tessd database was devdopei The database is 
structtffed ito korise ihe Itoiaiy data ffles, modeled files, derived measuranent, features, 

aod othar descriptive data. A master pot idenfificatioinmmbeF is used as fiie key to link these 



wo 03/088085 . PCT/USQ3/10665 

data elements with additional vessel data in existing databases. This structure permits additional 

databases to link to the query engine by adding a field with tlie master pSitffi la each record set, 

and developing an XML schema and DTD to link related data fields between the databases. It is 

scalable as a proof of concept, is consistent with Dublin Core cataloging structures, and requires 

minimal coding to provide access to additional databases - 

A visual query interface as described above was used to permit users to interact with the 
data using sketches or by selecting sample vessel shapes to augment text and metric search 
criteria to retrieve original and modeled data, and interactive 2D and 3D models. 

This project, successfidly implemented a powerfiil system of 3D modeling and analysis 
techniques to provide descriptive data and support research into vessel shape and structure. It is 
interesting to note that even measurements that are relatively coarse fiom a computer science 
modeling perspective offer significant improvements in accuraey fox cerainic researchers. The 
ability to search and cotDpare these accurate models of vessels offers new tools to cerainic vessel 
researchers.. 

Example 2 -LithicTocfi Manufacture cmdR^ttting . ■ : 

The gpal of the lithie refittifig pilot prcgect is the reconstruction of stbne. tool 
manufacturing activities within prehistoric archaeological sites and related cultural behaviors by 
identifying sequences of conjdiniBg lithie artifeGts. Refitting lifliic artifacts by manual trial and 
error is an accurate,, but highly labor-intensive, meHiod that requires the entire sanqile of artifacts 
to be present in a single lab. These conditions are not always possible given various antiquities 
restrictions. Automated .3D sur&ce matching of conjoinable artifacts will (frainatically enhance 
the efSciency of tUs valuable analytic mdhod and extmi His scbjpe of seainhes to collections 
&om different sites. Lithic refitting, or assembling fhe pieces of stone, broken apart by 
- pr^istorie people, has proven very usefiil in our atfempf to better understand the prehistoric 
•past The technique is especially usefiil in technological studies, taphonomic research, and 
spatial analysis. The 3DK project is developing software that "would reduce the time cost of 
refitting and facihtate development of an electronic database storing 3D images of lithic 
cofiectiotns available for study by researchers all over the worid. FIG. 46 describes the XML 
: sdie*na.Tised for iitiiies. 

Example S-3D Topography of Joint Swf aces 
The objective of the ''SD topography of j(a^ 
uaderstandingiof the abilities, Kmitations, and adaptations of our early humafl^^ceslors to make 
tools and walk upright by developing Monxediamcal mo^ls of miSttiipulative and locomotor 
behavior using 3D osteoloi^cal data. Use of calipers and visual kispeiGtion are inadequate to 
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capture tiie complex curvatures of 3D Joint surfaces and to control for body size differences in 

cross-species comparisons. This can be overcome by developing scalable'j"-<pjantitative models 

of reciprocal wrist, hand, and knee joint surfeces that will allow for comparative quantative 

analysis of the effect of surface area, curvature, and congruency on joint mechanics in extant and 

fossil apes and humans. 

Since the project inception, more than 600 bones representing the wrist and hand joints 

of humans, chimpanzees, gorillas, orangutans, and gibbons have been digitized to create a 

database that will eventually include approximately 1000 bones representing the vrast, hand, and 

knee joints of humans and apes. With an aim to better understand the fimctional moipholo^ of 

joiut suT&ees, the segmentation of features from a bone tbsi are of particular interest to a 

. physical anfluppologist such as joint surfaces may be automated, avoiding manual digitization of 

such features that is both time consuming and labor intensive. The surfece areas and curvatures 

of joint surfaces and the congniencies betweaa reciprocal joint sur&ces aire then quantified' and 

analyzed Static and dynamic models can be built i^ng digitized osteological data, tbgetho- wititi . 

musculoskeletal data from cadavers, and manipulative and positional bdiavioral data, to analyze 

the mechanics of manipulative and Ibcoinotor bdiavior. FIGi 47 shows the sdiema used for a 

bone. 

Featia-e Extraction from Volume Data 
Volume Segrnmtation Using WeibuUE-^ Fields 

According to another advantageous aspect of the invention, reg;ions. of volume data can 
be extracted for ejqjloiing the.inner structure of the volume data by perferming segmentation of 
voliune.data. Many tasks in volume visualization involye exploting the inner stnictures of 
. volume data. For example, a cell biolo^ majr be ioterested inihe stroctoie of the naic^t^bul^ 
spindle apparatus in an egg. The rapid increase k.jdata s^ skes repaired for coUectiBg images 
around the spindle ^paratu^ ^ well, as the j^or. si^ial to n£ftse ratio in &e data set, make it 
difficiilt to extract geometric features efS<aaatIy< . 

Segmentation of volinne data is a procegs sof voxel elassification that extracts regions by 
assigning the individual voxels to classes in such a way that these segmented regions posses the 
following properties: (1) voxels within the same region are homogeneous with respect to some 
characteristic (e.g.i gray value ortexfoie); and (2) voxds of nd^orii^ regions are significantly 
ifififerent with respast to the same chiffacteristie, 

The input data is on a 3D structured grid of verdces v(r, JyA:) , eadi assodafed with a 
scalar value. A voxel is coasidered as a k xie xk cube, ^and each voxel is assigned two values: 
oxp^tjaicy and standard deviation'(E-SD): . .'^^ index is used to estimate the 
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noise in a voxel, and to obtain more precise E-SD values for each voxel. The segmentation 

method has been tested using synthetic data as well as real volume d^fioni^£onfocal laser 

scanning microscope (CLSM). Analysis of this data shows distinct and defining regions in their 

E-SD plot. Under the guide of an E-SD plot, an object embedded in real arid simulated 3D data 

can be efficiently segmented. 

According to the volume segmciitalion method, the input data is on a 3D structured grid 

of vertices v{i,j,k), each associated with a scalar value, and a voxel is considered as a cube 

including k xk xk 3D structured points, called a k- voxel. Each 7c-voxel is assigned two values; 

expectancy and standard deviation (E-SD). The expectancy in a voxel relates to its mean/ and 

the standard deviation indicates the variability of the data within it. It is assumed that the E-SD 

values of voxels in a region are relatively homogeneous and different fi-om that in other regions. 

Many voxels have the same E-SD value. If one plots the frequency of voxels tht have the sanie 

E-SDi then some areas in E-SD domain will be dense and some spars©. This plot is called the E- 

SD field of the volume data. '.As will be apparent to those of skill in tfa& art, for a gjven volume 

data, the E-SD field depends on tbs size. of K^voxels selected, i.e., the value of k< 

A simple and efficient way to calculate the E-SD is to compute its average and tiie 

. sample standard deviation. However, noise malces it difficult to calculate the E-SD values 

accurately. Under this situation, the result of the E-SD plot is not stable and is dependent on a 

statistical model of the data. A number of statistical firameworks have been proposed to model 

image and volumetric data, "fhe ofasenred imajge pixels have been modeled ias Rayleigji 

distribution random variables with means depending on thai positioa A Gaussian-flmction was 

used for pixel relaxation labeling . Others instead proposed an ejqponential family of functions 

including Gaussian, Gamma, Rayldgh, and Possion to perform segmentation on 2D, 

According to the present method, a Weibull probabilistuj. fomeworfc ibr segmenfafion of 

Confoeal Laser Scanning Microscope (GS^SM) volume data is described. The Weibull 

xJistribution, fkst infrodnced in 1939 by Swedish engineer W. Weibull, builds on en^iiical 

grounds in statistical theory of the strength of materials. The WeibuH disMfailUoff includes' three . 

parameters, winch are described below. An- advantage of fiie Weibull distribution is that its 

kernel shape can be contioHed by^electing'di^fetieof paitoetiare for liie gray levels of fhe input 

volume data. ■ 

The following desciipEion discusses s^JatMy dtsttibutedofajects, tiie Weitull distribution 

sand the Weibull noise index, as weU as k>w to use &© W©S^ 

: algoritiims. 3D results -from coatrol data audtwo real GLSM wIUBie data sets afe presented. 
Spatialtjf Distributed Ob|eets. 

Ckjnsider volumetric data, r, is _ 

- ' r - m '^^.^^^^^ I '--: 5Bt /;^v--^--;-- : ----^^^^^^^ 
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given by v(/, k) or v(/>) , whose size is N=N^xNj,x N^. As used herein, distribution means 

the distribution of v over a certain interval {a,b]{a > 0) . The randonrvknable' J!fQ(v) is the 

number of points in a region QcF which have the value v, written as- A'q. The density or 

frequency /fj(v) of a random variable is defined as follows: 

where O,^^ = {(i, j, k)eQ\ v(i, J, k) = Vq} , and | Q | denotes die number of elements in Q. . 

Assume that a homogeneous segment can be mathematically specified usuig two criteria: 
(1) the relative constant of regional expectancy and (2) regional variance of the intensity. These 
criteria are expressed as follows: 

Definition 1. A. region Q. is called as a spatially distributed object (SDO), if the 
expectancy and standard deviation for each K-voxel A in Q are relatively constant, i.e. 

E[XJe(e„e,\andSD[XJe(d„d,), (2) 

where €,,62, if,, and cTj denote predefined constants, the random variable is defined 
as above. • ■ 

In general, the expressions of expectancy and standard deviation in a K-voxel are given 
as follows: 



EL^,] = r^,^y(w},^SD[XJ=l±- ^v\x,y,zy-E\XJ, (3) 

wiLere I A [denotes tiie nimiber of elanents in A. The SDOs are also call^ "agents" and regions 
of interest (ROI's). The goal of segmentation is to locate SDO's. The dioice of 6^ ,62 , rf„ and 
depmis on the E-SD field; . . 

However, if noise is present then equation (3) will not give accurate E-SD values. 
. CLSM. . data whichf , inherently includes noisc and has a poor signal to noise ratio resulting in 
. -these iiiaocui^ t^€5:,E-SD values. _ 

Confijcal LasCT Scanning Mieroscopy. (CSLSM) is a techaique &r obtaining high 
resQlutioii scans of optical slices througji a jfiiidE jqpecimen -Without having to cat the specimen 
-mechanically. Due to the precise lenses, -fllieJii^ coheicncy of the fllnminating laser beam, and 
.the confocal way of gathering bacfecattersd fi^ accocate fot^aang at ^edfic planar locatioaa 
can he achieved A typical optical section is hetween O.l^lQO [un . Scanning through the whole 
spedmen iieceJjy gives a fidl 31) projectitMi view of the specimen. Hiis technique! is very useful 
: not onfy fiie volumdric analysis of biolo^cal data?, but also because the 

; techm<pies used in *lstaiBajjg" .fliese spedmens (Le. l£f£ier excated dyes) increase the accuracy of 
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these images as compared to images obtained fi'om ordinary optical microscopes. Nevertheless, 

images are still noisy and blurred. Several sources of noise can be idenj^fiedJThgse include (i) 

thermal noise induced by the photomuJtiplier, (ii) photon-shot-noise, (iii) biological background 

(auto fluorescence), and (iv) unspecific staining. The quality of the image can be affected by a 

possible mismatch of refractive indices, tissue scattering, dye concentration inside the cell, and 

the histological methods used for staining. These factors contribute to a position-dependent 

noise and blurring, which makes the analysis of these images rather difficult. 

Statistical theory has been used for segmenting medical and biological data. This 
assumes that the data follows a distribution. Intensity values in a region have been assume^ to 
follow a Gaussian distribution. The Gaussian, Rayleigh, and Poisson distributions have been 
discussed separately. In our paper, before attranpting to segment volume data using statistical 
theory, we first analyze the distribution. FIG. 47B shows the distribution of a CLSM data (see 
FIG. 47A. The plot 600 shows the distributfon of the complete volume data (see FIG. 47A, 
wfaidi looks like the Poisson distribution. The plot 602 shows the distributioB at the brightest 
region in FIG. 47 A. This looks like a Gaussian distribution. The plot 604 illustrates the 
distribution in a 4-voxeL 

Weibull Distribution , 

Weibull distribution is defined as follows: 

■ where v> v,, a>0 is the shape parameter, b>0 is the scale parameter, and vq is the shift 
parameter (the minimum possible value of the random variable). In the CLSM data, the 
minimum possible density value is zao, i.e. v, = 0. Iherefore, it is assfumed that the shift 
parameter of the Weibull distribution v,, =0. HG. 48 gives the Weibull distribution (4) wilh 
.di£f^nt sh^e parameters a and scale parameter 6=1,2 and vo=0. The expectancy and the; 
de\iation of flie random vari^IeXobejdiig the WeibuH distribution ar^ 

£{Arj=6r(i+-li+v,,and fiD"[Jn=^>fr(l+^)-r(l+-)l, (5) 

« : L « <2j 

where the gamma function i5r(a) = J* t e''dt . It can be shovm that when.c =1.0, it is the 

Pdssion pd^ and when a=2.0, one has the Rayleigh pdf; and when a^S.O, it turns into the 
Gaussianpdf: When ^r»l, thedistribu£ioat€ndsto bea unifannpdf. Therefore, the WeibuU 
model is a suitable model to fit the histogram distribution of fliese volume data and the regions 
mtSmllkmi vfhose stB^^ controlled by 
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selecting a different a valiie. 

From equation (5), it is clear that the parameters a and b of .Ite-Wbitiull distribution 

depend on the expectancy and the standard deviation. We denote the ratio r = SDI E in equation 

(5) vdth V(, =0, and the relationship between r and the shape parameter a of its Weibull 

distribution is as follows: 

^2 ^ ^ ra+2^_ J 

or 

where the Beta function i&B{x,y)= \t (l-t)'^dt, a is the sh^e parameter of flie "WD and 

t = 1/a . From equation (6), it can be seen that fihe shape parameter of (he Weibull distribution 
is only dependent on ratio r in the E-SD plot. When t « 0.42 , ftxe liHS of eqiiation (^ reaches 

■ its maximum, which is near the value Q. 72. Unfortunately; the RHS of equation (6) is not a . 
monotoliic function of t (see FIG. 49). For each r, there are two roots, ia order to overcome ftus 
dfjfficulty, we first give some properties of the Weibull distributipii: 

Property 1 : For every j > 0 , the j-momeait of Weibull distributioa is: 
EX' ==b'r(l + s/a} . 

Property 2: If JC^,X2,---,X„sie independent distribution random variables, and follow 
the Weibull law, tfaestt 

Ly'xf ^ EX' for l"<J<oo,as 

» =w . .-. 

Weibull Noise Index 

' ' ' Removing noise or improving th.e sigoaf-to-noise ratio (S^]R) of a given iffl]age.;is m 
essentidl step ui segmmtation> especially in hig^ noise situations that can disriq>t (be sliape and 
lose the edge in£>iination ^lat defines the stnictate of objects. The traditional algorithms of 
denoising, saich as Gaussian fitter, reduce the noise but they do not mdntain the edge 
informatioiL When noise is removed, it is desirable, not only to smooth all of the fioraogenous 
regions tii^ contain noise, but also to ke^ theposition of boundaries, i.e., not to lose the edge ; 

■ tofoimalion that defines the structure of objects. 

Let Vj,V2,...,v., represent k- rniage points in a given K-voxel, It is assumed that the 
value of a voxel is characterized by the WeSttdl disti&tition. If equation (3) is used to calculate 



.55 
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the E-SD value, the results are not reliable due to noise, especially for a standard deviation. 

Therefore, one must find a way to distinguish whether or not the data disifibution in a /c-voexl is 

uniform. If it is not uniform, then what kind of noise is present? If few of the elements in a 

voxel are significantly larger and/or smaller than others, then these are called upper/lower noise. 

For example, in a 2-voxel, in which the set of the intensit)' at eight image points is {234, 52, 64, 

46, 50, 54, 62, 3}, the element 234 is much larger than others, and is called an upper noise. The 

element 3 is significantly less than others, and is called a lower noise. In order to classify the 

noise in a jc-voxel, an auxiliary function g(s) is introduced: 

= (7) 

where j e (-oo , oo ), v,. > O , l < /' < n , and « = k ' . By calculating equation (7) directly, we have 
the derivative of the function g(s) that satisfies: g'{s) ^0 for s >0 , > 0 for ^ < 0 , and 

k, k 

g'(s) = 0 ifandonlyif V, =Y2 = *" = v„. Also: g(0) = l, g-(<») =-l-/aad g-(--^) = — , where Arj 

n « 

is the number of the elements which are equal to the maximum, and k2 is the number of the 
elements which are equal to the minimum. It is obvious that l<k„k^£n,and ki+kj^n. 

(EX'^ t Bit t \ 

Using Property 2, we know that g{s) « ° ' 2' > where t^, -. s/a , and a is the 

Weibull distribution shape parameter in the K-voxet From the analysis above, the function 
t,Bit„t^)/2 reaches its unique positive maximum near 0.72 it i « 0.42. If Ai and h are small 
enough (^,,^2 <[0.14nJ^) and thereisa 5oSVUsh&atg(*o)==K).72,tbmw^ 

a^sJOAl (8) 

Definition 2: If > 0 sudi fhcat ^(j,) = 0.72 (if ^ > G.72« , set j, = oo)i-tien j, > 6 is 
called the Weibull uppw noise index. If jj>0 such that g(-iSj)=a72 (if Ai>0.72«, set 
=»-«), ftien >G is called the WdbuH lower noise index. In short are called tbe 
Weibull noise indices. 

These two parameters are used to determine the "goodness" of voxel distribution as 
follows (see FIG. 50):. 

For a K-voxel, if the Weibull upper noise index jj <L26 and the lower noise 
index ^2 > 1-26, tlien there is uppernoisem it. 

For a K-voxel, if the Weibtffi upper noise ittfex; j, >1J26 aa^ 



\x\ is (he maxinmnimteger'R^idi isless fii^nx 
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index 5^ < 1.26 , then fliCTe is lowra- Hoise in it. 

For a K-voxel, if the WeibuU upper noise index 5, <r.26- and the lower noise " 

index < 1.26 , then fhere is upper and lower noise in it. 

Segmentation Algorithm 

Based on the analysis above, the algorithm for volume data segmentation is as follows: 

Step 1: Given a k to detemiine the size of K-voxel, initialize the SDO's 
predefined constant in equation (2): e^> ei>0,d2>d^>0,. and the threshold of 
expectancy T^>6. 

Step 2: Consider the jffh K-voxel. Use bisection to- compute its Weibull noise 
index j,, and^^ which are the roots of the equation ^(j)=0.72, where g(s) is defined by 
equation (7). If fhoB is upper noise or Ibwer noise or both, tihen remove the noise 
directly (i,e. . delete the minimum or the maximum or both). Repeat at most jcic ^ J times 
to execute Step 2; 

Step 3 Calculate E-SD values; using equation (3). If the expectancy is larger than 
the threshold r,, add the K-voxel to list j4. If fhere are some jc-voxels which have not 
been dealt with, then go to Step 2. 

Step 4: Compute the frequency of the voxel in the list Ai- md create the E-SD 
plot. By usuig initial E-SD . values ^ >e, >0,tf2 kO, select tlie voxel which 
sa&s&es: e^>E>ei, d^'k.D'z.di, 

In this algorifhm, the flireshold is used for controlling the size of Iist.<4 above, and will 
cause the image to be rendered festear. The constant C in Stq) 2 is equal to 0.14.:. E-SD values 
(e^te^id^idi) aie obtained by moving a lectangle, called a window in the B-SD field, tbrou^ 
user interactiaB. 

The' algoridim desoibed above is simple and efEicienL Its average con^^ty is 
0(ZIogZ), wbereZ is thenHmba-of boxes, defined by X = 7//(k*) where// is the number of 
points in the volume data. As will be i^jpareut, if flie grains (or Kr-voxels) are coarser (i.c. k. is 
iarger)i fliaa the .atgoritibotn is more efiScieiit; however, the results of the segmentation iwU be 
coars«:also. With different k (ic =k k^^,k^ and are prespecified), deselected k: • 
best fits ih& given volume data; By this fitting approach, the number of SDOs in volume data is 
determined. The minimum ;^ which shows the number of SDOs in E-SD fiel4 wiU be chosen. 

Examples 

In fiHs section we iwill look at tvro ^exa^]pIes fflusfrating ffie pn^ 

■ -■■ . - • , .55. • . : . 
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segmentation. The first example examines artificial volume data generated usmg a WeibuU 
distributed random number with different parameters. The second exajspls^^ real CLSM 
data from two dififerent data sets, and demonstrates how this model can be used to segment such 
data. The hardware we used is a DeU Precision workstation 330, with P4 r,4GHz CPU, and 1- 
GBRAM. 

Controlled Experiment 

In Older to maike the e3q)eriment comparable, a controlled expeiiment is done in the 
following way: we first segment a noise-free volume data (see FIG. 52A) and. treat that 
segmentation as our reference, which includes a torus, an ellipsoid, and two defonned cube^ of 
different sizes but with the same shape parameter a. in a. 100x100x100 cube (see FIG. 52B). 
Every instance of a yolume data with added noise is then denoised using a Gaussian filter with cr 
= 1.3, and a WeibuU noise index with 2-voxels or 3-voxels, The targets are then segmented 
from these images and compared to the reference objects. The comparison, based on the 
segmented volume, is done by identifying the support function of the reference object and of tlie 
object segmented &6m a. volume data with added noise, denoted by and iS„, respectively. 
Th^ is iS^Cjc) -1 or 0, if A- is in the segmented objects or not, where % = {>-,»} . Then; the 
volnine deviation (error) of from in the vohmie data V is defined as follows: 



These deviations iEire calculated to produce numbep that are coirparable acarpss different 
noise levels. Several levels of noise have been added to the test volvme data to show the 
robustness of the filter. The noise that is added to every image point is a Weibull. distribution 
.jWth different scal^ parain^eis b: y==jwbi^55,(:^bli^-JOf'h '^'^ Rpidom wagUtble X i^ thfi 
uniform distribution, a is the -Weibull sh^e parameter and C is a constant fer each object 
Finally, the volume ercor from the simulated volume data is plotted in FIG. 5L As depicted in 
FIG. 51, the voluma error conespondii^ to the /WeibuU noise index is signifieantly lower 
conqjared to those that result from applying a Gaussian filter. Segmentation results are shown in 
FIG. 52. 

Figure 6(i) shows the Weibull E-SD fidd of the noise-added volmne data^ w^ 
parameter b=lQ.O . Below we define where the ooloi^ in m E-SD field coirespond to 
fiiequency. We denote fay Nie^ ib& nvanjb^ of flie xr-yoxels whose expeotaajcy K e and standard 
deviation is (f. Tlxe color at point (e, tO in. the E-SD field is detamaiiKd lQ^ log(iV(€,<f)) - We set 
the color at point (e, d) as follows: 
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'i?G5(255,255,255), 0 <log(A^(e,J)) <0.5 

RGB(25S,0fi), M.S<\og{N{e,d)<l.O 

J?GS(0,255,0), 1 .0 < \og{N{e,d)) < 1 -5 

color{e,d)^\RGB{0,0,255), 1.5<log(iV(e,if))<2.0 
J?G5(255,0,255), 2.0 < \os{N{e,d) < 2.5 

RGB(q,255,255}, 2.5 <lQg(N(e,dy)<3.0 
i?Gff (255,255,0), 3.0 <. log{N(e,dy) 

The colors in the E-SD field below have the same meanings. Next, the rectangle, called a 
window (see FIG. 6(1)), has two small circles at its left-top and right-bottom vertex, which give 
the values (e^ ,e^;d2,d^) to define the range of expectancy and standard deviation of a SDO. v 

Figure 6(c) shows a slice of the noise volume data with fiie Weibull scale parameter 
6=5.0. The segmentation results, using our method and threshold method with Gaussian filter, 
are given- in Figure 6(d) and 6(e). AUhoiigh it has lost some detail information, such as the 
. defeimiation in the two cubes compared with the reference Figure 6(b), the segmentation using 
. our method keeps the number' and shape of objects. In contrast, segmentation performed by 
. using threshold methods with Gaussian filter lost the number and shape of the objects.! 

Figure 6(f). shows,^a. slice of the noise, volume data with the Weibull scale parameter 
6=10.0, and the segmentation results from this noise data using our method and tlie threshold 
method witli a Gaussian filter are given in Figure 6(g) and 6(h). Figure 6(i) shows the Weibull 
E-SD field of this noise volume data, and illustrates that the three different components present. 
Using our segmentation method maintains the number and sh^e of the objects, but using 
direshold techniques with Gr^ussian filter Mis in segmenting the objects. 
CLSMData 

In this example we use data from two different CLSM data sets designed for 
investigating the ioeiotic spindle wifiiin a mouse egg. The eggs were viewed on the Leiea TCS 
. NT confbcal microscope. Multq)le las^ allow Sai simultaneous imaging oftbe DAPI ( Argon 
UV [363 am]) ALBXA S6& Kryptoti[568TmJ) ttaoropkoit-ls^^ U&ng a 63x 

w^otgecfiwe, images were scaaned l:fetweertfl.2rtl.4 jun mtervals 'al<Hig the z-axiis, 0.154 jmi 
slcm^i^X^&ids atid collected Atbi^ flie volume of fiie CTitire egg. 

Hgure 7 shows a CaiSM test bed for oar WeibuirE^D modeling scheme firam diff^ 
experiments. 3te' data in Figure 7(a) is coUected usmg a Kryptoa laser and bigjdighte regions 
.targeted witii an antibody to anti-ct-tubulin at tfie upper left and brighter regions through the egg. 

; The image in Figure 7(a) shows a meiotic met^hase U arrested egg, and is composed of 124 
sKces within a 512 x 5 12 matrix contains a gray-ievel fiom 0 to 255. The size of the data shown 

: & Figure 7(b) is 512x512x146, «ind has the same gray-level range as in Figure 7(a). In Figure 
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7(b) the egg has a meiotic spindle at the upper left and the remains of the primary polar body at 

the bottom right. 

Figure 8 shows the E-SD plots of Figme 9. The colors have the same meanings as in 

Figure 6(i). hi Figure 8(a), the size of the window is (170, 237; 4, 27) corresponding to the 

segmentation of the data shown in Figure 9(a). hi Figure 8(b), the size of the window is (183, 

255; 2, 30) corresponding to the segmentation of the data shown in Figure 9(b). We set the 

:^threshoId at -=34, and r=2. The segmentation in Figure 9 (a) includes 3902 voxels, 

• and 1 1308 voxels in Figure 9(b). 

Figure 7. Test bed for our Weibull E-SD modeling scheme. The images come fronvfwo 

different data sets. In botti images there is only one cell and its spindle is at the up part of the 

cell (labeled by biologist). The size of image (a) is ; , its gray-level is Som 0 to 

255. The size of image (b) is 5i2xsi2x 146, it has the same gray-level as (a), 

... Figure 8. Hie Wdbull B-Sp fields fcom the data shown in Figure 7. The coldrs have the 

same mean as in Figure 6(i). In 8(a), the size of window is (178, 237; 4,.27} corresponding to the 

segmentation shown in Figure 9(a). Thethreshold - is 34, which creates a blank 

on the left side. There are two clear regions. In 8(b), the size of ivindow is (167,255; 2, 30) 

. conespoudiQg to the segoientation shown in Figure 9(b). The threshold ■ i s 34, 

which creates a blank on the left side. There are two clear regions too. . 

Figure 9. The spindle segmentation eorresponding to Figure 7. The segnentation in 9(a) 

includes 3902 boxes, and 11308 boxes in 9(fa), 

What gives rise to such clear regions m Wdbull E-SD field shown, iii Figiire 8 is 

Tinknown. A pfeuable ©cplanatioii is that each region corresponds to a SDO. When we move the 

wmdow on Figmre 8(a). to the small region, an area of cell perimeter in FiguoBe 7(a) is segmented 

^(see Figure, 10). . .- . -• . .-. . ' 

gonclusion . 

, We have proposed a coarse-grain approach for the . segmentation of an object -from 

volume data based on .tiie We&uH E-SD field. . We have shown that one can decide the noise . 

index in a K-voxel by nsiag Weibnll law^aad use the E-SD fiejd as a guide to segmient an 

object We have cionsistejitIy>demonstrate4 this, spproaxiti on' coalrolled as well as on real 

volimie data. . ■ • . 

. A new ^proach for segmenting 3D TOxd. data s^^ 
automafio and .joonsists of two phases.? ;3^ns^ tl>e initii^ se^watalion is accda^Kdted fay:voxel 
labelmg using statisttical majdmnm .KkefilKXjd cstimaficm ted^q^^ of the present 

.^58 ..... ■ ■ ^ :■ 
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approach is that the type probability distribution function is not required a priori. A multi- 
parameter distribution that includes a variety of widely applicable distributions is utilized. The - 
second phase refines the segmentation by the use of a new greedy connected component labeling 
(GCCL). The overall effectiveness of this new approach is illustrated with data examples of 
multichan-nel laser scanning confocal microscope (LSCM) images where the structure of GFAP 
(a protein) and nuclei, and the geometry of an electrode implant are extracted. 

Multichannel laser scarmiiig confocal microscopy (LSCM) equipped with a specific 
laser/filter combination to collect multiple fluorescence signals in a single scan is applied 
increasingly in experimental biological investigation. In this application, the confocal 3D 
images are collected firom two channels. The stack of 2D images (each 2D image referred to as 
a "slice") make up the volumetric data. A typical section is from 10~IOOnim. It is composed of 
different sh^ed structures, for example, tree-shaped GFAP (a protein expressed by astrocyte 
. cells), and blob-shaped nuclei (labeled with a DAPI stain) and sheet-shsped electrode implants. 
The- segmentation of these regions can yield important biological infbnnation. For example, 
segmentation and subsegment determination of GFAP volume allows for quantitative analysis of 
the immune, response to the. implanted electrode. In LSCM in partieuliar, measurement of the 
relative positions of regions labeled with different cells/implants can provide an insig^ht into their 
inter-fimctional relationship. Due to image noise and different sh^es af flie cells/implants to be 
segpi^ed, eonsiderable e£fort is required to develop acGurate image segmentation for localizing 
the structures. 

Recently, several i»obabilistic firameworks to determine algorithms for the segm^tation 
in an imi^e have been given. It is now well known that to obtmn effident atgorilhms, the 
statistical properties of voxels should be taken into account Menibets of the exponential fainily 
includes Gaussiai^ Gamma, Rayleigh, Poisson, and other finailiar distributions that have been 
used to model real data. Indeed, there arc several kinds of 2D images for -vdiich the pixel values 
are ^correetly described by such statistical law. Gaussian distribution is widely used to 
characteriz© the probabilistic feature of random variables of an image, ome ultrasonic medical 
images liave been modeled using Rayleigh's law. Astronomical images have been presented 
using Poisson distflbutioir and compared to the classical linear intercorrelation technique. 
Previously, it has been assumed tiiat the probability density fimction ^dj) of flie random fields, 
wMch describe tire gray levels in the input ima^e; belongs to flie expotieritial femily. However, 
:a j3f^^^vdfe a single distribution sh^e will 1^^ the use of tbe segjoaeutation qsproach, and asks 
users to distitigBish the difife^atf j)dfi of input linages, fc; aiccordance witfi fte'^esent , we use 
itbe W^ifaull j?^^ whose kernel shape can be con^ 

WIiCTl a statistical model-based meHjod is used, both voxel and its nd^ibois ^uld be 
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considered, to estimate the parameters of the pdf for each voxel. Then, a histogram, called the 
Weibull a-b histogram, is generated. Using voxel labelmg guided by the JV^ihuHa-i histogram, 
the initial segmentation result can be obtained. However, due to the unavoidable noise in LSCM 
images, the initial segmentation is quite coarse, e.g. this results in many isolated small regions. 
In order to overcome the problem, we introduce a new algorithm of connected component 
labeling (CCL), called Greedy Connected Conqjonent Labeling (GCCL), to delete the unwanted 
small regions. CCL algorithms segment the domain of a binaiy image into partition that 
corresponds to connected components. CCL is one of the fundamental segmentation techniques 
. used in image processing. The classical CCL algorithms pay more attention to the cos^ of 
memory but not time complexity. The most common approach to avoid making several passes 
over flie image is to process the image in a sequential order and keep an equivalence table that 
stores the information about the equivalence of labels as-signed to the voxels fijom the saine 
component. Others have used two passes in which the image is scaimed in raster order, and the 
size of equivalence table is bounded to 0(N) for a Nx N2D raster image. Due to the increase in 
tiie size of a 3D image, GCCL takes into account the computation time cost as well ^as the 
smaller size of the equivalence table. In the GCCL algorithm according to the present invention, 
one pass for a 3D image is used, the size of local equivalence table is Ofihe widtk of connected 
component), and tlie time cost is O^^ with low constant, where iVis the number of voxels in the 
image. 

The algorithm according to the present invention has been tested using LSCM volume 
data shown in FIG. 60. Figure GRl(a) shows a LSCM image of GFAP-labeled astrocytes and 
DAPI stained nuclei^ in which the green indicates regions targeted with GFAP-labeled 
astrocytes, and the red regions identify DAPI, which consistently tai-gets nuclei cells. Figure 
1(b) diows a GFAP and implant LSCM image, in which the red targets electrode implant 
devices and the green GFAP, FIG, 60A and 60B are rend^ed using maximum intensity 
projection (MIP). 

Statistical Image Model 

We will now in&oduce Weifaull statistical modeling and describe how to create a Weibull 
a-b histogram and how to ob^ an initial se^nentation. 
' Weibull Modeling 
Consider a 3D image Joigarazed as a sl^dc of 2D transverse dicesi 
is d(aioted by/, whose siz© is %:xi^xiVjW3ieiie the intoigity of eadi voxel ^ e lis gjvea 
vfi;/; k), Ths WabuIlprobabiHly density fimction of gray level of eacli voxel is j^ven: 

by 
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where v > vo; a > 0; 6 > 0; vo > 0, and v is the gray level of the voxel, a is-a shape parameter, b is 



a scale parameter, and vq is a shift parameter (flie minimum possible value of the random 
variable). This triparametric distribution was introduced in 1939 by Swedish engineer, W 
WeibuU, on empirical grounds in the statistical theory of the strength of materials. It can be 
shown that when a = 1 .0, it is the Possion pdf; when a = 2.0, one has the Rayleigji pd£; and when 
a = 3.0, it turns to the Gaussian pdf. When a » 1, the distribution tends to a uniform pdt 
Therefore, the WeibuU model is a suitable model to fit the histogram distribution of ihese 
images and the regions in them whose statistical properties are unknown since fhe kernel shape 
can be controlled by selecting a different a value. FIG. 61 shows tiie Weifaull disttibution 
equation (1) wifli different shape parameters ^i, the scale parameter 6 = 1.2 and the shift vo = 0. 
For LSCM imaging, the minimum possible intensity value is zero, i.e. vo = Q. Therefore, in the 
discussion that foUws, we assume that the shift parameter of the WeibuU distribution satisfies vo 
= 0. 

Assume that the observed image / is composed of two zones: the target having One (or 
several) simple connected region(s) and fhe background. Under this assimiption, the target ia 
the image is completely defined by a binary image T- {T (i, j, k) \ (i, j, k) e 1} that defines a 
certain shape for the target so that T (i, J, k) is equal to one within the target and to zero 
elsewhere^ Thus, the target in the image is the region: Q.T={(i.j, k) e I\ T(i,j. k) = l}. 

The purpose of segmentation is, therefore, to esiEimate the binary image T&ii th©. target in 
the image. Without more a priori knowledge about the taiget, the Tnaximinri likelihood 
estimation is first introduced.to confute the Weibjill parameters. 

MgYtmiim T.ilreHhnod Rstitnatinn 

The parafneter estimates obtamied using the maxttnom likelihood tecbmque ar@ unique, 
aiid as the of fhe san^fe increases^ the esthnates statistically .t^oach. the tme values of the 
san^e. Let vl; -v2; ^ 0 ; vn repress n voxels, of image< It is assumed ^t the in-tensity of 
eac& voxel is diaracterized by the Wdbuli distiibutioa The lilselihood function associated with 
the sample is the joint den-sity of n random variables, and flins is a fimction of tiie unknown 
WeibiUl distribution parameters (a; b). Hie. likdihood. fimctioti for a sample under these 
assimiptions is given by the expression 



The parameter estimates are detennined by taking fiie partial derivafives of the logarithm 
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of the likelihood function with respect to a and b and equating the resulting expressions to zero. 

The system of equations obtained by differentiating the log likeUhood-feac.t^n.fpr a sample is 

given by 

In order to find the parameter a from equations (3), we first in- troduce a function g(a) 
defined by equation (3). For any sample vl;v2; jifjzSf; ;vn with vi > 0(1 • i • n) , we can directly 
compute that tlie derivative that satisfies g 0 (a) >0, and so g(a) is a monotonic as- cending 
function with shape parameter a, and lim a ! 0 + g(a) = ( ¥, and lima ! Ppg(a) = kln(max(vi).) V 1 n 
llhi=l In(vi) > 0 , where k , I is the number of the elements, which reach die maximum in the 
set f vl;v2; ^0 ivn g . Therefore, for any sample vl;v2; (z!^(Zi ;vn with vi > 0(1 -i -n) , g(a) has 
one and only one positive root. Figure 3 shows the plot of fiinction g(a) for different samples. 
Once a is determihed, this value is inserted into equation (3) and b is cal- eulated directly. From 
, equation (3), it is easy fo see that b is the a-moment of sample vl;v2; ;vn with vi > 0(1 • i • 
n) , and a indicates the deviation of this sample. The less deviation, the larger the root of 
equation g(a) =0 (see FIG. 61). By the analysis above, we can use the bisection algoritbifn to 
solve equations (3) and get the Weibull parameters. 
Vbxel labeling by Weibull a-b Histogram 

Once die WeibuU model is obtained, the segmentation problem amounts to assigning 
labels to each voxel in the volume. A straight- forward way is to label voxels as the target or the 
background by maximizing the individual IScelihood fitiK;tioiL This appioach is called ML 
classiiG^, which is equivalent to a multiple thresholdiiig method. Usually, this metiiod may not 
achieve good perfbnnance since there is a lack of local nei^orhood infonnatioii to be used to 
mate a good decision. Thecefore, we incorporate the local neighboitood infonnation at a voxel 
into a labeling piocedare, thus improving the segmaitation perfbnnance. Suppose v i; j;k is the 
mkodty of a voxet 0; j^) 2 1 in tbe image I,, w i; j;k is a aze d £ d £ dvmsdow centred at (i; 
j;10 foa: maximum likelQiood estimation, where d is an odd integer ^a£@r tiian 1. We regard w i; 
j;k as the local region while we calculate &e Weibull parametccs fac tiie voxel (i] j;k) using 
equation (3). 

"When Ihe intensity value ranges fi-om 0 to 255, then the value of tiie WeibtiU scale 
p^^ameterb is also fi^om 0 to 255 by equation (3). On other hand, fhe more umfofitn the region 
surrounding a voxel is, the larger the WeibuU shape parameter a becomes. Bxperimentally,. we 
set the upper boundaiy of the Weibull shape parameter a to be lOQ. The isize of the window has 
hsSneDce on tiie calculation of the WeibuH paramefets. Hie wJodow should be large enough fo 
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allow enough local mf ormation to be involved m the computation of the Weibull parametersTor 
the voxel. Fuitliermore, using a larger window in the computation ofidje^pa^iGgters increases 
the smoothing effect However, smoothing the local area might hide some abrupt changes in the 
local region. Also, a large window may require sig- nificant processing time. Weighing the pros 
and cons, we choose a3x3x3 or 5x5x5 window for estimating the Weibull parameters and 
have experimental data to back that choice. 

A classical histogram is a statistical gr^h counting the frequency of occuirence of each 
gray level in the image or in part of an image. We extend this idea and define a histogram in the 
Weibull a-b domain. First, the Weibull shape parameter a and scale parameter b for each yoxel 
are calculated. Second, for each Weibull parameter pair (a;b) 2 [0;100] £ [0;255J, count the 
number of voxels with this parameter pair. Here two issu« arise in computing the frequency for 
each parameter pair. One is that the step of the Weibull a-b donmn is (1,1). The other is that we 
set a low boundary for the scale parameter h, since b is the a-moment of the intensity sample 
sunrounding a voxel We should identify the target, with higher intensity, not the background 
with lower intensity. Last, v^e have the fitiquency for each parameter pair logarithmized, and. 
plotted against &at pair, and colored as the legend shown in FIG. 62. The Weibull. a-b histogram 
gives us a global description of the distribution of the uniform regions across intensity levels. 
We use a movable rectangle in the histogram to locate the range of the colored peak zone by 
moving its top-left or bottom-right vertex, and select all the voxels with die Weibull a-b 
parameter histogram being in this rectangle. 

An advantage of segmentation using the Weibull a-b parameter histogram, is that local 
information and global information are both taken into account in determining the segmentation, 
whereas in traditional histogram approaches only global information is con- sidered. 
Segmentation Refinement Using GCCL 

Although LSCM images have a much high accuracy, these images are still noisy and 
blwed. Several sources of noise can be identified: fliennal noise induced by the phototiiultiplieE, 
photon-sbot-noise, biological background (autofluoiescence), and imsjpedfic sfaitung. Hie 
qudify of the image d^mds also on possible mismatdi of refiactive mdiees, oa tissue 
scattmng, on dye concentration inside the cell, and &e histological methods used fiw staining 
methods. These factors contribute to a position-dq)endeht noise ^^^m^ 
segmeaatation results using voxel labelling by fte Weibull a^b histogram can be quite coarse and 
may lead to the probleni that there are many isolated small segmented ccahponehts, as shovm in 
.FI0. 64A. This is due to two reasons. First, a thresholding based segmentation is a binary 
representatioii of a region that either includes or excludes a gjveii voxel as beang in or out of the 
region. Second, the voxel UOielliag cannot distfaguish flie fcffgets fiom the noise. On the other 
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hand, voxel labelling can not show the sh^e of ceUs and implants segmented and the 
relationship among them, because voxel labelling only uses flie in-formatton-amimd a voxel. 

We wish to correct these problems by deleting the unwanted small regions and finding 
structural and quantitative information in an image using connected compoiient labelling (CCL), 
which seg-ment the domain of a binary image into partitions that correspond to connected 
components. Here, two voxels are 6-connected if at most one of their 3D coordinates differs by 
- 1, 18-connected if at most two coordinates differ by 1, and 26-cQn0ected if all three coordinates 
are allowed to differ. A 6/18/26-connected path througji the 3D image is a sequence of 6/1 8/26- 
connected voxels. A 6/18/26-connected component in a binary image is a subset in whid» for 
every two voxels there is a 6/18/26-comiected path .between them. The partitiordog is 
represented by an 3D image in which all voxels that lie in the saine connected component have 
the same voxel value. Distinct voxel values are assigned to distinctly connectied component. It is 
clear that only the target voxels are affected by this labelling i.e. the background voxels remain 
unchanged. 

Hierarchy Frame For a Connected Component . . 

For a connected component G, a hierarchy fiame H C (r) rooted fyym voxel r is defined 
as apartitioning of C into hierarchy f hl(r); li2(r); hn(r) g satisfying 1. hl(r)= f r g ; 2. All 

voxels adjacent to voxels in hierarchy hi (r); (z = 2; ^ ; j 1) are in hierarchies hi \ l(r);hi (r) 
and Ai+l(r); 3. All voxels adjacent to voxels in hierarchy An(r) ^e in hierar-chies hn \ l(r) and 
. An(r). where n is the total number of hierarchies,and is simply the depth of C. maxf § hi (r) jg is 
called the width of C, where j ^ j denotes the number of the elements in a set Figure 5 shows a 
2D 8-connected component which is labeled by the set f A; B;C; D; E; F; G; H; I; J; K; L; M; N; 
O J P; 2 g • If the root is A, then the hierarchy Same for this connected cort^dnent is hi (A) ~iA 
g , h2(A) =f B; ^ g , h3(A) = f C; F; Q;^ g , h4(A) = f D; Bg, h 5(A) = f H; Jg., h.6{K) = f K; 
L; M g , and h7CA) = f P; O; Ng '. Therefore, flie deptti of G is 7, and the ^dth is 4. 

<3GCL • ■ ^ ^ . 

For a bifi^ imagfe T^ ifi^kt GCCL :&ds"^- unJafiefed vdxet callisd ^ root i, it does not . 
stop uatii aH voxels in JFFcoimected to r titougji a (5/f8/26-connected jpafifi are labeled. The 
: prooediireoflabelling a connected eoiiE^ohentC by GGCL^ 
roQteclfixHfir: piir impleraeirfatidii jo^ as adynasiic 

rfist -OOL), aMidifferent com-ponmts are assigned to ififf^^ 

GCGL is ffaat GCCL may rq)eat to label a voxel in C. In order to avoid the e3q)ensive 
Qpraation of DL, sudi as adding a unique node to a DL, we use 4-vaIued flags to determine the 
. value of a voxel in T. Let t(i; j; k) dmote the value at voxel (v, jj fc) 2 T, and set t(i; j; k) to be 
dfii^ 0 or 1 or 2 or 3. If t(i; j; k>= 0, thm the voxel j; k) beloi^ to the background and is not 
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changed. When t(i; k) = 1, the voxel (i; j; k) belongs to the target, it is not la-beled, and can be 

a root to form a new component or a new node to be added to the DE-Whaa t(i; j; k) = 2, the 

voxel (i; y, k) can be a new node to be added to the.DL, but can not be a root to fonn a new 

component When t(i; j; k) = 3, it means that the voxel (i; j; k) has been added to a DL, neither as 

a new root nor as a new node. The GCCL algorithm is as follows: 

GGCL Algorithm: 

Read, Zand set t(i,j,k) = I, if (i,J.k) ear; ' 
while t(i, y, it) = = 1 loop 
. /* the first loop*/ 

Generate a DL, denoted by to store the coimected component whose root is /c), 

and a temporaiy DL, denoted by tdl, to store middle results. Two counters, denoted 6y 

numl and nuni2, represent the increase otdl and tdl, respectively, and are iiiitialized to 

2ero. Add voxel ^) into rdl; 

while tdl 9t 0, and voxel (/ ,m, n) e tdl loop 

/* flie second loop*/ 

i£f(/ .OT, 7z) = = 1 or 2 then 

while all the 6/18/26-connected voxels o{(l-,m, n), denoted by (7± 1, m± 
1, n ± 1), and /(/± 1, m ± 1, K ± 1) = Hoop 
/*the third loop*/ 

nvtml++; . . " 

Add voxel (/± 1, m ± 1, « d: 1) into fjfl; 
Setr(/±I,m±l,H±l) = 2; 
End loop 
. Sett(l,m, «) = 3; 
Add. OT. «) into £fl- 
num2++; 

. End if 

Remove voxel ,m, n) from the M; 
if huml ^ num2 then 
Reset the f£fl; . 

End loop 
End loop 

Here, four issues arise in the GCCL aigoriflmi. One is that for a comiected component 
Cits root is the voxel (lO; jO; kO), where kO = mm (i; j;k>,2 C k, jO = min (i; j;kO) 2 C j, and iO = 
. min (i; jO;kO) 2 C i. The second is that we use three primitive qperstions of DL: 

Add(3?): It creates a new node for contaitung a voxeL IlJfinode is then pushed on the DL 
..ofvoxel x at 0(1) time cost - -. 

: Remove(x)^ This operation fiist.moves the.pobtM' of DLcto -the root of DLi^ 
finds the voxel X m fheDL; If xis m'DL,lhm delete, If not, -go 1^ In genaral, if 

&ei:e are « nodes in fhe DL, tins operation costs Q(n). Inifais algorithni, the voxel to be removed 
is always the root of the DL, therefore, this operation has 0(I)^eost. 
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Connected Component 




Data 


Data Size 


Amount Largest - 


— Time(s") 


GFAP & Nuclei 


29.9MB 


2514 87536 


4 


GFAP & Implant 


63.9MB 


3115 234976 


8 



Table 4: The efficiency of the GCCL algorithm Data Data Size 
Connected Coiiq)onent Tiine(s) Amount Largest GFAP & Nuclei 
29.9MB 2514 87536 4 GFAP & Implant 63.9MB 3115 234976 8 

. Reset: It reinitializes the DL and returns the root of the DL in 0(1) cost. The third issue 
is that the number of times the &st loop is excuted m the GCCL algoritlim is the number of the 
connected components in the image. The number of ttnaes in the second loop is dependent on 'the 
number of voxels in a connected component, and the third is a constant either 6, 18 or 26. 
Therefore, the time complexity of the GCCL algorithm is 0(N) with small constant, where A'' is 
the number of the voxels m the hnage. The last issue is that the size of local equivalence table is 
0(the width of connected component), because each local equivalence table is the DL tdl in . 
GCCL algo-rithm. hi FIG. 64, the right image is the voxel labelling result by means of Weibull 
S'b histogram, includes 2514 26-eonnected com-ponents, and the left one is the resuUs of the 
GCCL. Table 4 shows GCCL on diflferent data. This was run on PIV 1.7GHz with 4GB RAM 
Window 2000, As shown m FIGs. 65A and 65B, the resulting surface from our segmentation 
caii be very coarse. This is due to the noise in the images and thiesholdingbased segmentation 
methods. We smooth connect components using convolution with Weibull kernel whose 
parameters are determined by the equation (3). FIG. 65 shows the smoothing results of 
connected components. 

FIG. 60A shows LSCM volume data composed of 58 slices with_an in-plaue 512 x 512 
pixels, and IGO x 100 x 57:8 mm field of view. Afiercomputetion of the; Weibull shb histogram 
^wn in FIG. 62, we move flie top-ieft vetle^ to (24, 0) anja &e bdftdin^ii^ to (255, 100). 
After the initial segmentation using voxel labelling by the Weibiill a-* Ms-'togi^m, the 
corresponding result TO is given in HG. 64A. Usmg GCCL fer TO an^ setting tiie threshold of 
the number of voxels in a connected coiap<me3A 200, we can get the sdgifltaifeifion resiilte shown 
in FIQ. 66. FIG. 67 shows the efifect of the top-left vertex on Weibull a-£> histogram on the 
connection, for glia cells in LSCM image FIG. 60B. By adjusting: the position of VCTtex on 
Weibull a-6 histogram, the connections between adjacent astrocytes can be more clearly defined, 
allowing for the potential idaitification of insfividual astrocyte. FIG. 6t)A is obtained, when the 
toprlefi; vertex is moved to (3 1,0), and FIG. 60B to (24,0). r 

, As pother example, we consider a GFAP and implant LSCM .vdbmp data shown in 
FIQ. 64B, FIG. 64B is 127 slices with an in-plane 512 x 512 itiatrix. FIG. 68 shows the 
segmentation results of ttie GFAP and inqilant data.^^e have tested bur algorithm vidfk 69 >. 
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different LSCM volume data. As illustrated here, we have observed that the overall approach 

can be very efifective at segmenting LSCM images. 

In this paper, we have presented a new statistical modeling of vol-inne data to segment a 

target of interest, and a GCCL algorithm to refine the initial segmentation from the Weibull 

statistical model-ing. This method, as illustrated by pilot application iri LSCM im-ages analysis, 

is Citable of segmenting the structures within data and can be applied to real problems such as 

those encountered in tissue segmentation. One of the remaining limitations of the present 

approach is that it is still semi-automatic and consequently requires the intervention and 

expertise of the user. It would be desirable to move in die di-rection of a more fiilly automatic 

segmentation procedure. 

Conclusion 

The above-described system and method of the present invention possesses numerous 
advantages as described herein and in the. referenced appendices. Additional advantages and 
modifications wiU readily occur to those skilled in the art. Therefore, the invention in its 
broader aspects is not limited to the specific details, representative devices, and illustrative 
examples shown and described. Accordingly, departures may be made from such details 
without departing from the spirit or scope of the general inventive concept. 

Appends I 

. <?xm| yersion="1.0" ??r 
- <xmlpot> 
:i <pot> 

<ld>6</id> . . 
<pot^Id>130S2</potJd> 
<lln6ari.unlts>centlmeters</Mnear_unJts> 
<volume^units>c^tfnietei9.</vQlunfi^units> 
<areaL.units>cent!m4ters</are^_unlts> 
• <referenoe_:^stem> . 
<xl>l;0</xl> 

■ •• <;x2>0;O-</}j2>: ^ • ■ ■ ^-'w.-- A.- .•- ■ ■ ■ - 

<x3>0.0</x3> 
<x4>0.O</x4> 
;. .-. <x5>0.O</xS> 
<x6>1.0</x6> 
<x7>0.0</x7> 

<x8>0.O<yx8> ^ . ' 

. <x9>0.0</x9> 

<xio>o.o</xio> 
<xll>1.0</xll> 
<xl2>0.0</xl2> 
<xl3>0.0</xl3> 
<X14>0.«</Xl4> 
<xl5>0.0</xa5> 



</reference_system > 
^ <measured_value> ^ 

<compactness>0.858936</compactness> 

<vesseLheight>15.2221</vessel_height> — 

<min_diameter>12.0058</min_diameter> 

<max_diameter>16.9216</max_diameter> 

<surface_area>788,948</surface_area> 

<volume>2344.09</volume> 

<symmetry> 1.0</symmetry> 
</measu red_va I ue > 
^ <form> 

< vesseLtype > Restricted </vessel_type > 
<contour_type>Inflected</contour_tyj3e> 

</form> 
^ < representation > 

< ra w_data > n u II </raw_data > 
<model_data>nun</modeLdata> 

<raw_url>vidya.prism.asu.edu/disk2/KDI/Data/pot;s/iatO& 
2.ply</raw_url> 
<modeLurI>NA</model_url> 

< thum bna il_u rl > http : //vidya.prism.asu.edu/kdi/potth umbn 

ail/302S2.gif<:/thumbnail_url> 
<thumbnail_data>nuH</thumbnail_data> 
</representation> 
j: <profile_cun/e> . ; 

<curvejd> l</curv.ej.d> 
<pronie_reference> 

<pronie_pQs}tion>Arfaftrary</proflle_position> 

<a>1.0</a> 

<b>0.0</b> 

<G>0.O</C> 

<d>0.O</d> 
</profile_reference> 
- <representation> 

< raw_data > null </faw_data> 
<model_data>nult</modeLdata> 
<raw_url>null</raw_url> 
<model_uri>null</model_urt> . 
</repFesentation>< : . /i r- ,.v 

i <2one> : ..V \ . V ??- -i 

- <rlfn> 

<diameter>12a79»4</dfameter> ^ 
</nm> •. . I 

.;^<neGk>.., •••<vr ^--^-^i. .. 

<neck_diameter>0.0</neck_diameters- ' > 
<h€fght>O.0</height> 
<corner_pointl>0.0</comer_polntl> 
<comer_point2>0.0</cx)mer_polnt2> ^ 
</neck> 

- <orifice> 

<diameter>0.0</drameter> 
</orifice> 
^ <body> 
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<ciiameter>O.0</diameter> 

<height>0.0</height> 

<start_pointi > 0.0 </start_pointl> 

<enci_pointl>0.O</end_pointl> 

<start_point2>0.0</start_point2> 

<6nd_point2>0.0</end_point2> 

</body> 
ji <base> 

<sta rt_point> 0.0 </start_point> 

<end_point>0.0</end_point> 

<curvature>0.0</curvature> 

</baSe> 
</zone> 
— <feature_point> 

<point_id>l</point_id> 

<point_type> EP</point_type> 

< t_pa ra m eter> 0 . 0 </t _parameter> 
</feature_point> 
<feature_point> 

<point_id>2</point_id> 

<point_type>EP</point_type> 

<t_parameter>55.0</t_parameter> 

</feature_point> 

<feature_point> 
<pointjd>3</polntjd> 
<point_tvpe>EP</point_type> 

< t_pa rameter>110.0 </t_parameter> 
</feature_point> 
<feature_polnt> 

<polnt_id>4</pointJd > 

< pQint_type> CP </poi nt_type> 
<t_parameter>3.0</t_parameter?- 

</feature_point> 
± <feature_point> 

<point_id>5</point_id> 

< poin t_type > CP </poi n t_type> 
<t_parameter> 107.0 </t_parameter> 

</feature^o(nt> . 

<feature_po(nt> 
<point_id>6</point_id> 
<po|nt_type>iP</point_type> 
<t_parameter>-l.O</t_parameter> 

</feature _point> 
i <feature_point> 

< poi nt_id > 7 </poin t_id > 
<point_type>lP </pointJ:ype> 
<t_parameter>10.0</t_paFameter> 

</feature_polnt> . 
± <feature_point> : . . - 

, <point_jd>8</pointjd> 
<point_type>IP</point.type> 
<t_j3arameter>102.0</t_paraiT»eter> 
</feature_point> 
^ <feature_point> 

<po'nitJd>9</pointJd> 
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<point_type>IP</point_type> 
<t_parameter>108.0</t_parameter> 
</feature_point> 

< / p ro fi I e_cu rve > 
</pot> 
<potdata> 

<Spedmen>13082</Specimen> 

< Project> RPMS </Project> 
<Form>Jar</Form> 
<Site>U:8:24</Site> 
<Type>Plain Smudged </Type> 
<Phase>Gila</Phase> 
<BeginDate>1320</BeginDate> 
<6nc!:pate>l450</EndDate> 

; : ■.^<^F6ATTYPE>Z7</FEArrYP£> 
. <StRATUM > 238?</STRATUM> 

■ <STRATA>0</STRATA> 
<SVOL>0.0</SVOL> 
<SUBFEATU>238.04</SUBf=S\TU> 
<SUBTYPE>6</SUBPi'PE> ' • . ' ^ ' ■ 

<SqREEN > 0 </SCREEN> 
. <ARTIFACT>:l</ARTIFACT> 
<S1CSTR^A>.3</S1CSTREA> . 
<slcstre.alabel>pJain/smudged</slcstEealabel> ■ 
<VESSPART>2</VESSPARTi» , 
<vesspartlabel>Neck -r Jar</vesspartlabel> 
<LUSTER> 0 </LUSTER> 
<speccount>l</speccount> 
</potdata?- 
:/xmipot> 
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Abstract 

A new approach for segmenting 3D voxd data, sets is presented. It 
is semt-automatic and consists of two phases. Firet; (he initial seg- 
mentation is accomplished by voxd. Isdbdiiig using statistical max- 
imum lilcelihodd estimation techniques. The novelty of the present 
approach is that the type probahiiity disttibution fimctian is not re- 
quired a priori. A multi-parameter distribadon wiiich includes a 
The second 



variety of widely applicable distcibutioiis is ufilizei 
phase refines the segmentation by the use of a new 
neoted cotnponetit labeling (GCCL). Hie overall eJ 
this hew approach is illusbated mdi data examples of multichan- 
nel laser scanning confocal microscope (LSC&Q images where die 
straotnre of GFAP (a protein) and nuclei, and the geometry of an 
electrode implant are extracted. 

CR Categories: t4.6 [Image Processing and Computer Ver- 
sion-]: segmentation — Region growing, partitioning; M.IO (Tro- 
• age- Processing and Computer Version]: Image Representation — 
Statistical, Volumetric 
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Figure 1: Test bed for our- volumetric segmentation scheme, the 
right one, denoted by (a), shows an ESCM image of GFAP-labeled 
asfrocjftes and DAPI stained nuclei^ in which tile green highlights 
regions targeted with GPAP labeled astrocytes, and the red regions 
identify DAPi; whiA consistently targets nuclei cells. The left one, 
denoted by(b), shows a tSCM magt of GPAP-Iabeled astrocytes 
end decttode implant; in ^duch the led tacgets implant and the 
gteon GFAP.Both ate iBn<fctBd by *HFin 3D. 



1 Intraduction 



laser scamung confocal microscopy (tSCM) 
equipped mths specific lasec'fllter oombiaation to collect multiple 
fluoiescence -signals m & sisgle scan is applied increasmgly. in ex- 
perimental biological investigation. la our^plicatioi^ the confocal 
. 3D nnages collected from two channelj. The stacfc of 2D images 
- (eaA 2D image lefened to as a ^ic^ make up die votumetric da^ 
ptazdan et aL 200tJ. A typical section is fiomlO ~ 100 fim. It is 
rfofdiflerentj" " ' - - - 



itDCtet (labeled wi& a DAPI stain) and tbtet-d^ed electrode ■ 
plaats. Thesegmtfotatitxi of Aesen^onscanyieldinipcntantU^ 
logical fafoimation fSacA et sL 2000]. For example segipcntatioo 

tatiw aaabtsts^f the umnme Mspmtse to.aeiDipltmled ekctode 
.pyitnaretat. 1999]. rnLSCMfnpatBeiga^ me asm ' Bii e ii tof tte 
iel^ve positions of legtons labded with diffeicnt edib/lnqilants 



can provide an insight mto their inter functional relationship. Due 
to image noise and different shapes of the cells/'iraplants to be seg- 
mented, considerable effort is reqnhed to develop accurate image 
segmentation for localizing the structures. 

Recently, several probabilistic frameworks to determine algo- 
rithms for the segmentation in an image have been given (Chanda 
et aL 1988][Lee 1998]Ea>e8naud 1999]iCh^arabotty and Duncan 
19991jjA et aL 2001}. K is now well k&Eywn that to obtain efBcient 
algori^ims, the sta&deal properties of voxels should foe taken mtc 
. account. Mambas of die exponendal ftmily includes. Gaussian, 
Qamma, R^ei^Poissoo. and odier fehuliar distribntioas [Chcs- 
naiidi999Jtfaat faavebeeni]sedtoraodef real data. Indeed, there 
w sewtd loads of 2D images for which the pad values aiB cor- 
iccfly descdbed by sudi statisdcal lacK. Gausaan distiibafion fLi 
etal. 2001] is rnddy used to characterize the probabilistic feahve 
of random variables of ah image. In [Chanda et al. 1988], some 
Tilhasonic medical images are modelled osing Rayleigh's law. In 
fLee 1993], Astronomical images have been presented using Pois- 
son disiribiitipn and coinpat^ to theclassical linear intercorielation 
tBcbmqpe. Si [Oiesnaud 1999], the authors assumed that the prob- 
,9hSBty 4soa^fm<Am fields, which describe 

i)ie gfs^ ievds in the inpiit image, hidings to Qie exponential 
itr. HowE^ api^rtritb a single distttoition sfc^e p:,i et al. 2001] 
rnfflliinit.tiie usd bf beseg^ and asks nseis to 

distiagui^ die difTeteiil^s of iiptt imi^es. In fiiis paper, wc use 
the WeibnUpi^ whose kHiid sbajie can be controlled by selecting 



od is used, both voxel and 
its aei^ibois diould be esnsidefed, to estimate the paiametas of 
" ^ dfaesiiaoa i999t|pjetaL 2001]. llieo,a 
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histogtam, called the WeibuIla-6 histogram, is generated Udag 
voxel labeling guided by the WeibnlIa-6 histogram, the initial seg- 
mentation result can be obtained. However^ due to the unavoidable 
noise in LSCM images [Manders et al. 1 992], the initial segmenta- 
tion is quite coarse, e.g. this results in many Isolated small regions. 
In order to overcome tiie prdblem, we introduce a new algorithm 
of connected component labeling (CCL), called Greedy Gomiected 
Component Labeling (C5CCL), to delete the unwanted small re- 
gions. CCL. algorithms segment tfie domain of a binary image into . 
partition that corresponds to connected components. CCL is one of 
&a Amdamentat se;graentation tedmiques used in im^e processing 
[tifoga and Gabbouj 19971. The clas^cal CCL algorithms pay more 
attentioil to &e cost of memory but not toe con^Iexlty. The most 

■ comiaonapproadi to avoid makingseYeral passes over the i^iage is 
to process the image in a sequential order and ke^ an equivalence 
able that stores the hifoimatlon about flic equivalence of labels as- 
signed to tiie voxels fcam the same cothponent: In (Khanna 200 1], 
two passes hi which the image is sicanned in raster order are used, 
and the size of equivalence table is bounded toO(JV) fora/irxA:2D 
raster image. Due to the mcrease in (he size of 3D image, GOCL 
takes into account the computation tim;: cost as wellt'as the smaHer 

■ size of the equJvalenee table. In our GCCL algorithm^ one pass 
for a 3b image is used^ the size of local equivalence table i!0(the 
v4M ofcomectedeompoitmi, and flus time costisO(W) wife low 
constat^ wbecei^ is fte number (if voxels ia ifee iinage. 

We have toted our aijifoiithm udng LSCM voiunie data shown 
in I^e 1, ti@m lCa>sho«« iUSCM fmags of QFAP-labded 

' a^oK^ and DAPI staiaed iincla, }n wWA Ifae gteen in<acBtes 
tegJoiis targeted vri(h\GPAP-Iafaered astroigptes, and fte red regions 
identify DAPI, which consistently targets nuoiri cells. Kgore 1(b) 
stayrs a GFAF and impltet LSCM bilise, in which fee red tsi^ 

.^^dosfitp&SapIart 

al20dlL .. V - 

In &sdon It we introliice Welbidlia^stidal iefaHid&S andi^" 
SGnbie fasw to create a Wetbulla-A histogrim and fa<»r to ob^ 
initial se^aitatioa. la Section HI, we discuss s^mentation tw-' 
.-&einsattising .GjGGL. This is fin^er.extai^ to tialKnislnife the 
icoiSaiyiBdiiqirbdUdba^ 

wei^l;^jfee;alg9iifem tercatljSCSyC Images. F&ally, fiitsie wtxk 
Bfecuwetf'inSfeationV. . ■• ' 



. 2 Staltistical Image Model ,, . . 
2.1 Wabull ModdHng 

Oinsid*s3D&as^foiliaftizedasast^ 
and the gfM pomtse* iffdawted Ig^ vi*iaeeaze «5 it 
-wfaewllteintensitydf <ach 
IheWdbuffpmbabffity density ltocfii»yMi)o^^ 



Figure J: Plot of {unctiang(a) defined by equation (3) for four 
samples which are generated by Weibull pdf will* = 1 .2, and vq = 
0,0 anda = 0.7,2,3. or 10. respectively. And 0.7. 2.1. 3.05 and 9.8 
are the roots of g{a) = 0 corresponding to fee samples. 



voxelis given by tSuhir 1997i[Hirose and Hideo 1996] 

'<')=f(?)-'-K^)T 



paramster (the minimum possible value of the random variable). 
This tiipirametrie distribution was introduced m 1939 by Swedish 
enginCCT, W Weibull, on empirical grounds iii the statisa'ca] theory 
of fee strengfe of materials. It can he shown that whera = 1.0, it 
is fee fosslonpdf; and when a = 2.0, one has. fee Rayleigh;>4^ 
and whffli« = 3.6, it tums to fee Qs!»smD.pdf, When.a:i» I, fee 
distribution tends to a vai&sBa^. Therefore^ the Weibull model 
is a suitable model to fit the histogram distribution of feese usages 
and feese rq^ons in feem whose sta^cal ptopetdes are unknown 
Since the knnel shape can be controlled by selecting di|ferehb 
value. Kgure 2 ghres Wdbull distribuflon equation (I) wife differ- 
ent the shape parameteFse. and fee scale parameter6= 1.2 and fee 
sUftvb =::0. For LSCM imaging^ fee: minimum possible faitcnsity 
value is zero, i.e. v^j =0. "Iliefefore, in fee nest part: of feis pa- 
per, "we assume feat fee shift parameter of fee Weibull distribution 
sat)8fiesM> = 0. . 

We assume tliat fee observed image/ is composed of two zones: 
fee target having one (or several) smqite; (^niiected re^on(s) anid 
feebadcground Under feis assumption, fee t»get in fee imige is 
completely defined a bmary iiriage 7= {r{i, /,i)l(i,/,<r) Si} 
feat defines a certami shape fiirfee ta^so thk te equal to 

one \rifein.fee target and to zero dsewh^ I%us, fee taiget in fee 



iniage>feej<(«iw:jQr.^ {(U*J 
. .Tjhepdip(»e«&A^;^eata)iQnt^.ja><7e^^ 
iiw«ftrfertte'.l|iigeS&% ?W^dBtifeoie-a^qri knovd- 
■edtae aiwat'^"'''^'^^^ ' 
esuD»£ii)a) 



•2.2 ■ Maximutn Likelihood Estimation} . 



tecirfque are unique^ and as fee size of fee sample, jnereases, the 

estimate'statisticaUy z^prqadi tfae1ni9 vaft»^ 0f^ Let 
. •'h'^i' reptesentn vakz of jinai^ "ft is^^ssnmaifeat fefe in-. 

teod;^ of eadi voxelis laiaiactsrized^by feeWatM'^stribarion. 

Ihe Mdihood functfon assooated wfe the kui^t 

sitv o£a random variables^ and thus is a fimction of fee unknown: 
. \^>idistijj^on-pa!B^^ . :; 

' * ' 'OTIff.ti<)&'isglV&fa^fe^^ 
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Figure 4: The Weibull parametera-6 histograni for flie LSCM im- 
age ^o\;«m in Figure l(a)«i<i the color legend, v/her^(a,b) is the 
number of voxels whidi satisfy ^ye^bull parameters are fzand b 



The parametei estimates are detennined bx taking die partial deriva- 
tives of the logarithm of the likelihood function with respect tm 
Slid £ and equating the resulting expressions to zero. The system 
of e<BUUions obtained by dilftienliating tiie log likelihood fiincdofli 
for a sample is given by [Hiiose and Hideo 1 996] 

»<">*^^-:5'^-'-:-».'-[:(j:f)|"'- 

■ • . (3> 

- litoiderto fiiid'thepaiametBra fix>mi»jiiattons (3), wefirstiii* 
troduee a ianction£(<i) defined by equa^ (3). For any sample 
VI , »^ , - - • » v« vntb: v» > 0(1 < i < »} , we can directly compute tfiat 
<iiBdenvativetliatsatisfiB6g'(a} >0,andsog(a) isamonotonicas- 
cending function with shape.parameteio, and linv-.(i+^a) = 
and limfl-4-.g(a) =fchi(max(vr)) - ^17=1 In(v() > 0 , where fc> 1 
is the number of ^e elements, which reach the maximum in Hie 
set {vi,V2,---,v„}. Therefore, for any samplevi,v2, --,v„ witti 
V( > 0(1 <i^n), g(a} has one and only one positive root. Fig- 
ure 3- shows the plot of functioqg(a) for difTerent samples. Once a 
is determined, this value is inserted into equation (3) ancfc is cal- 
culated directly. From equation (3), it is easy to see ttiafe is the 
a-moment of sample vi,i^,- with > 0(1 < t , and a 
indicates the deviation of fliis sample. The less deviation,.^ laigcr 
tiieTootofequation^a) bO (see Figure 7). .By &e analysts tbave, 
we can use the bisectian dgoti&unto solve eqaaBona (3) and get the 
V^bibnU pdcametet£ 

.2.3 Vopcet.laiieriHig tiy VVeibuir a-b Histcferam 

Once tiK. WdMl modd is obtamed, &i st^mintation probletD 
- aau>imts to aissfiiiiii^.Iabels to eac^ voxd m flie V9lnm6. A stnaght^ ; 
■ . j^wtDDd, w^ is p fiibd voxbls as the target or tfc teidt;gnfund bj 
. 'iqax±n^mg;&ei imimn liliiModd fusistioa. tfiEs'^ioa&b & 
ifSlaiML cM^|q; wiu(% is a^^vderit fa airndt^ flsdiiolifing , 
mediod^ Usuatlji ms mdSiiod j^tuAa^^c goad perfimnaiioe' 
since fliere is a lade of locd n^AoAobd infcsmation to be nsed- 
to make a-£iE>od iikcisi<ni{Wao^^aL 2001]. .TteiBfore,'we inccir ". 
,.porateae.to<?4ndpiboitt^ . 

"v^fe As intaasi^ of a^^^ €/-ti fte&agcXiJ^iFilsA" 

es&natUon, vAeretf is an odd mfeeeE-greaiBr fiiaa 1. We icgisd 

ibrtiw TO*d(i,/,it} iisingx:qaaJion (3). 

. 'Wdieii einteiisi^WiE laages ftorn 0 to 5ti55, ftca^e.vdtK! of 



• I— • 
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example of hierarchy ftame for a 2D image 



enough local information to be involved ia the computation of the 
Weibidl parameters for ttie voxel.. Furthemiore, using a larger win- 
. -dowin the cpmputatioa of the patanieters increases the smoothing 
. effect However, smoothing the local area might hide some abrupt 
changes in the local-region. Also, a.Iaig^ window may require sig- 
nificant processhig time. Weighing Ifad pros and cons, we choose a 
3 x3x 3 or5x5 x 5 window for estimating the WeibuIJ parameters 
and have ffitperimental data to back that choice. 

A-'classical histogram is. a statistical graph counting the fr^ 
queni:^ ofoccoirence of each gra;^ levd in the image or in part of 
an image [Ritter and Wilson 2000]. We extend this idea and define 
a histogram in the Weibiilla-i' domain. First, the Weibull shape 
parameter o and scale parameter i for each voxel aie calculated 
Second, for each Weibull parameter pair(a, 6) e [0, 100] x [0,2S5j, 
coupt lie number. of voxels with this parameter pair. Here two is- 
sues arbe in cenqiuting Qie frequency tov each parameter pair. One 
is ^t the step of the Weibu11a-£, domain is (l,!>. Ttie other is 
UiatwBset.a low boimdarjr for. the scale parametsiA, siuceiis the 
armomoit of ^e intensity, sample surrounding a voxel. We should 
idendl^ fte taiget with higiier intoisity, not the baclcgrouad with 
. lower intensi^. Last, we have the frequency for each parameter 
pairlogaiifiunized, and plotted against that pair, and colored as die 
legead'Shown in Figo^ 4. The Weibulla-i> histogram, gives us a 
, global description of ^e (gstribntion of the uniform regions across 
ineaisity levels. .We use a movable Te^:Cangle in the histogram to 
locat? tte range of flie colored p.eak zone by movmg its tc^eft die 
hottom-ri^t vertex, and s^ect-all tiie voxels with &e WeibuIb-& 
pa^metcrjiiisto^ani being in this.rectangle. 

An iut^Bt^e ofs^entadon die Weibulla-4 pafiane>' 
ter Ustc^at^ fiiot JociEft iafatio^oh and global infbmiaiioa are ■ 
bodi taked^ii^ accoimt li debimiinu^ the segmentation, whereas. ' 
-in tsiditiottal iii$]iG;gFam i^ptqac^ cwdy global infoimatioii is con- 



3 r^^^meni^tf^ Refinement Using GCCL 

AlflMugh LSCliC Iniages have4i miicli high accuracy, these images 
aresGUtttrisjraiidldmicl Sev^ 

Otfeiiiidis^ iL 1992] : difaibial fioise indnced by the photomnitt 
~ pli^ pbo(aitdiot-4»is(^l)td[o^eal bidtgioimd (autoSuorescenceX 
and unspec^ stainiiig.' Hier'^ialify of file image depends also on . 
possible mismatsh of refractive indices, on tissue scattering, on dye . 
co^c^dii^d^itiMde tiiecdl^ methods used for 

; jstaiifingiaisai«is> .Thae fe^ts cohtabute to a position-d^endent 
rTtoise Siidttu^nng. Iterefor^ segmentation results using voxel la- 
.bd&^by a»Wabtflla-* histogram can be qms coarse and may 
lead to fte iwoaem that feere are many isolated small segmented • 
. oal^<Mie^ 3:.^wn aI^&e€{a): This is dne to two reasons. 
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Fiist, a thre holdini+a <; ^mentation is a binaiy representatio 
of a region tha either mcludss or occludes a given t oxel as being 
in or out of the teg on Se.,oii(l, the voxel labelling can not disCin 
guish the taigete trom the noise On the other hand, voxel labelling 
ean not how the shap of cells and implants segmented and the 
relationship among them, because voxel labelling only uses the in- 
formation around a voxel. 
We wish to correct these probtems by deleting ths unwanted 
■ small regions and finding structural and quantitative information in 
an image osin^ connected component labelling (CCL), wiiicli seg- 
ment the domain of a binaiy image into partitions that correspond to 
connected components. Here, two voxels are 6-connected if at most 



coordinates (U&f by 1 . and 26-connected if ait three ciioniinates at 
allowed to difiec A£/lS/26-coiinected path through the 3D image 
is a sequence of mS/26-oooaeeted wxeU. A 6/I8/2fr«onnecled 
component in a binary image is a subset in wtach for eveiy two 
' voxels &ere is a 6/1 S/26-aaaiiected path between them. Tlie parti- 
tioning is represented by an 31) image in wfaieh all voxels that lie in 
the same eoimect^ cotc^oneitt have the same voxel value. Distinct 
voxel values '^ assigned to distinctly connected component It is 
clear that onfy die target voxels are affected by diislabdliugi.e. the 



34 H^rard^ .Frame For' a Connected Contponent 



fianie£(;(r) molBd 
6om voxel r is defined as a pattitbnuw of C hOo tueraidiy 



• X M vwads atga^t to voxels in hieiiautiy %(r), (r = 
2, — ,B - i) are in liiaarcliiesft/.li (rj, i,-(r) and Ai+f {rj; 



where « is file total number of iiraiacdues^and is .^inqily ibey 
4cp& '0f C ,mmipii(t)\} is caned the widtJi . of :G, where 
it^ f aeaotBS Ac number of ftc daneiife in a ««. : ligtsie 5 



:._X^*t^^^t4^tJ>3: , ........ . . . 

..:tiiehieiah% &taie&r ffiis cotmected cdBmm^ttAiL&i 




GFAP labeled astracyte (a) before smoothing ,aiid (b) af 
ing, and a nuclei (c) befbie smoodiing and (d) after smo 



3.2 GCCL 

Forabinaiy imager, when GCGt finds a unlabeled voxel, called a 
mot t, k does not Oae unfil all vexsls ii^; connected to r though 
a 6/1 8/26-coHnefited pafii are labeled. The procedure of labelling 
a connected oomponentC by GCCL is to find the hierarchy ftame 
Hq{t) rooted fromr. In our implementation of GCCL, we regard 
a connected component as a dynariiie list (DL), and different com- 
ponents are assigned to diffaent DLs. However, the main problem 
with GCCL is that GCCL may repeat tci label a voxel inC. In order 
to avoid the expensive operation of DL^. such as adding a unique 
node to a DL, we use 4-V8lued flags to determine the value of & 
voxel in T. Ut t{i,j,k) denote the value.at voxel (/,/,*) e T, and 
set t{i,J,k) to be ddier 0 or I or 2 or 3. m{ij,k) = 0, then the 
vnKt&y,^} bdongs totitebadtgroohdaad is not changed. When 
MhjA = t, the voxel (l,J,k) belongs to the target it is not la^ 



TnodetqteaddedtofteD^tetcaaiiilteawottofe. , 

conipoiient Wiwuft/.iij ==3i UTO&Mm JKMi(l,J,K^ has. 
beeii»td«lta aDI^neifiieras anenr roots0F.as aaeviiode. The 
GCXX algoiiiiim is as fidioTO: 

GCCLAl^ritfain; 

/* the first loop*/ / • ; : : 

<i«iraite a p£^ denoted bi?^«^ to store ffi^H^ 
ponent \iiiose n)ot is{», J, Aj, aiid^^^ 
.; . . byti^ tb store iniddte results. Two raiotttte, denoted by - 
inimi andmian^ r^tesent thfriii««ase ofeT and td!, te- 
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Table I: the efficiency of the GCCL algorithm 



Data 


Data Size 


Connected Component ^ , 
Amowit Largrat' ' ~™c(s> 


QFAF&Nudei 




2514 . .87536 .. 4' 


OFAF jjs^iaut 




Sits. Z^m 8 



/* the second loop*/ 

irf(/,m,n)==lor2then 

whUe all the 6/1 8/26-connected wxds ofi[/,»>,ii}, 
denoted by (f ± l,m± l.u ± 1), and «(/±.l,in± 
l,n±l) = lloop 
/*tbe third loop*/ 

Add voxel (/±.l,;«±l,/i± l).inl»/<ft 
Set<(/±I,<jj±l,B±l>!=2; 
• End loop 

Seti(/,m,iJ«3-; ' 
Add (/,m,n) into <//; ... 
. . nuinJ-H^ 
: Endif 

Reiiiove voxel (OVi) fiwn tfie 
if nuials4inHn2UieB , / . 
RewtlheifciS; .' 



The fiiini issue is Oat fiie number of times the &st loop is excuted 
in the GCCL algorithin isdienumberofthe connected cotnponents 
in the image; The numberof times in the second loop is dependent 
•on the number of vomls in a connected component, and the third 
is a constant either 6, 18 or 26. Tlierefore, the time^ complexly of 
the GCCL algorithm is Oljf} with small constant, wfieteAf is the. 
number of the voxels in tfie image, the last issue is that the size 
of local equivalence table isO{the width of comaded component, 
because each local equivalence table is the UUdl in GCCL algo- 
rithm. In Figure 6, the right image is Oie voxel labelling result by 
ineatis afWeibvUtf^hisiogim^ meludes 2S.1426^nnected 
IMneats,anddie'le&onet5 fte results of tfa« OCCX. "Dbte 1 shows 
OCOL on difierent data. This was lua on ptV U7GHz witii 4QB 
RAM Wimtow 2000; 
As sfNwn inFigw« 7{a) and 7{c>,tiMiesal&ig«iiiaib£^^$om our 
— '-'-•-j can be very coarse. TU( is due to Sid^ioise in the 



otedfoa a»^igd& We anooth 
n witb.Wetibiiil Iseaiel whose 



Hersi fourissues arise in the GCCL algoritfun. One is liiat £ir 
a connected componaitC, its root is the y<s)sA(l<>Jo,kt), M*eie 

The secoiid fc that we use three primitive opcradom 

• Add{x^ ft cieatBS a ne5(( node OTntafaiing a voxeL Th6 
- i;aodeistteiipushedontlieDLofvoj«katO(I)tinw 

• "R^ovcfzy This operation fiistiiwvcs deira 
ft© locrf of DI, iMtiiat finds flie Voicel X 
.teDL; tiien dd«e, Knotgo-drod^.,^^ 
tfiei?$«sff nodes in Oe DL, das ope^ 

. V . ilgori^ tiie vrntcl to b« lenioved is afw^'flu coot of the - 
I>I^ftaefer^tfu$operafioiilii8sO(t)~cost - 



4 Results 

We have aitwjfy described ttie theory beJiind Weibnll statistical 
awiddHie and GCGI» 

eieexan^efcaSetofi^qSi^witumedatashoTO 
Bffiie 1(a) is 6bii^oscd<)f58sircesw^ 
e^ aiaiOOx ig0 k 57.8 /zM jdd of Afl&- computation of 
. Ife Wohollii^ hiStograni shoi*n fa Rgute 4, we move ihS ia^ 
Mstfex top4;^ to (255, tOO>, Aiier the 

s^meaitaftm nsing voxel 1*^!% 1^ ft© Wtibulb-i his- 
&^«^ te.<>«B5qp«»fiaf tesultTfl^^ 

<XXX.fx Ip and setting Hie thresKold of the ftumber of voxds in 
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-FigunB 10: The segmenfaiioa results of the GFAP and implant 
LSCM image Figure l(b>, (a) two channel data, and(b) die OFAF- 
laibeled astrocytes, and (c) die electrode implant segmented 



a connected component 200, we can get the segmentation results 
shown in Figure 8. Figure 9 shows the effect of the top-left vertex 
on Weibu'll a-b histogram on the connection for glia cells in LSCM 
image Figure 1(b). By adjustinE the position of vertex on WeibuU 
a-b histogram, the connections between adjacent astrocytes can be 

more clearly de&ied, allowing for the potential identification of m- 
dividual astrocyte. Figure 9(a) is obtained, when die top-left vertex 

is moved to (3 1,0), and Figure 9(b) to (24,0). 
As another exainple, we consider a GFAP and implant LSCM 

volume data shown iii EigiKe 1(b). FieiOs 1(b) is 127 slices wi& 

an in-plane ^12 X 512 matrbc. Figiae 10 shows die segmeqtatioa 

results of die€iFAP and implant'£& .. 

We have tested our algoiitjinr widi 69, diffisrenf.LSCM volume 

data. As illustrated here, we IiaVe b&loVed thU tfii^overaU ap 

CM be vtijr effective at s«?|fneitting ESCM imagra, - 

5 Conclusion 

In diis pape^-^^liaxfeffVeSented- a jiw i^aSs&tikiiio^^ugi^voV? ■ • 
mne data to ie^ei^ k tatlgei of i%r$s^ tm^lGC^i^^ . T 
refine die initfai^ segmaat^n fibm ^« We^ii^ i^a^s^ iiodetT;; ~ 
' ing. Ttus metiKi4 as iUustratedby: pilot.appti^dn ia tiSCM im- 
ages analysis, is capable of segmqitkg Ihf stimc^es 
and can be applied to real probiedui sodt^is ffiow epcq^ht^oi in - . 
tissue s^meat^tioo. ■ - ■ ■ :^ ■ ■ ■ 

One of the letiuUiwg linutations'df the prraeat it>e<oaicii &''ttn^ .v 
it is sdllsen^-automatic and oi>d^eqaeiifly iicgtdrEs ffie utfaii^ttidS 
and expertise of the user. It -would be de^ble tosMve iiHtxe dl- 
lecfioit of a mors tally automatic seigniec t a t i o n procedoFe^ 
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segmentation will be coarser also. Therefore, flie size of tbe^ -fcyVOxeL. needs 
experimentation to find the optimum defined by the user (cell biologists m our case). 

3. Examples 

In this section we look at two examples illustrating the proposed m«thod for 
segmentation. The first example uses an artificial volume, data generated using the 
Weibull distributed random number with different parameters. The second example uses 
real LCM data from two experiments, and demonstrates how this model can be used to 
segment such data. The hardware we used is a Dell Precision workstation 330, with P4 
l,4GHz CPU, and l-GB ROM. 
3.1. Cond'olted Experiment 

3.1.1 SiiButafing WeibuU Distributed Etandom Nambers 
We start with a randoin number, :e, which c&infes flom M identical distribution (in tide 




figure 5. RaMomly g«iaating.WeibuIl distdbution Using (7) 
with the scale parameter b=l^ three dktiiict:$bi^paramei^ a==0.7S,a=s3 and a=iO. 
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range from 0 to 1). We apply die transformation equation (9) belew~t&-get- a new 
independent random number which has a Weibul! distribution with a mean and variance 
that depends on the values of a and 6 [15]. Figure 5 shows random generation of a 
WeibuII distribution with the scale parameter b=1.2 and three distinct shape parameters 
a=0.75,a=3 and 0=10. 

y = [-bln(l-x)Y\ (9), 
3.1J2. Generating 3D WeibuII Distributed Volume Data 

In order to control our experiment and make it easy to render flie volume data, we 
generated two concentric spheres in a cube of size 100x100x100. The center of the 
sphere is (50,25,25). The radii of the external and flie internal spheres are 30 and 10, 
respectively. Tfie generated volume data has three parts, which all follow the WeibuII 
distribution. Their shape parametCTS are 0.75, 3.0 and 10.0 for the cube, the external and 
ttie internal sphere and their scale parameters are 1.2, 16.58 and 63.55, respectively. 
- Figure 6 shows the 3D WeibuII distributed volume data. 
3.13 Control Experiment Results 

Figure 7 shows flie E-SD plot of tibe 3D control volume data described above. First, 
we explain what the colors in an E-SD plot mean. We denote N(e^ the numW of the k- 
voxeb whose expedtahcy is e and standard deviation is J. The color at point (e, d) in E- 
SD plane is determined by the logarithm ofN(e/£), i.e. U)g(N(e,d}} . We set the eolor at 
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point (e, d) in E-SD plane as follows: ■ "-"^--^ • 

' SGBi255,255,255), 0 < log(N(e,d)) < 0.5 
RGB(25SA0). 05 < \og(N(e,d) < 1 .0 
J?G8(0.255.0), 1 .0 < log(N{e,d}) < 1 .5 

RGB(Qfi,255), 1.5 <\ag{N{e,d)) <2.(i (10) 
*G5(255A255), 2.Q<]oe(N(e,cl)<2.5 
i?GB(0,255^55), 2.5 S log(Ar{c,rf)) < 3.0 
RGB(255,255,0X 3.0^1og(A^(e,rf)) 

The colors in tiie E-SD plot bdow have the same meanings. Second, the rectangle, called ^ 
a window (see Figure 7), has two small circles at its left-top and right-bottom vertex, 
which give the values (ej,e,;cf,,ar,) to define the expectancy and standard deviation of the 
SDO. 

It is easy to find that the cone part on the left in Figure 7 corresponds to the part of 
the cube, which does not contain the spheres, the disc at the middle to the external sphere^ 
and the small region on flie ri^t-top to the internal sphere. Figure 7 shows that it is 
impossible to distinguish the outside sphere from the cube, if we just use mean field 
theory [10], because their projections on the E axis are overlapped. 

Figure 8 gives the segmentation of the contrpl 3P Weibull distributed volume data. 
Figure 8(a) shows the cube. Figure 8(b) ^ves flie external sphere, and Figure 8(c) 
illustrates the segmentation of the internal sphwe. The colpr in Figure.8 is RGB(155,v, OX . 
wjiere V is the value at tiie point (x, J/, z). 

Due to the tuitnber of ^-voxels mcreasiog exponentially cwresponcfiog to the 
color, we call the separate area in E-SD plot a finger. By equation (6), we can conclude 
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that the voxels from spheres follow the Weifaull distribution, in -"ijdiich-tha. shape 
parameters are relatively constant. Similar results can be derived for the-^xtemal sphere 

and internal oae. 




(b) the partial volume data with two co-cenfiic sphoies 
Figure 6. Generatittg 3D WeibuU disUibuted volume data m&i distioet shape and scale parameter 
Id. (a), tfaeblue part (cube) with ^.75 and b=12, thip red part (outside sphere) with a=3 .0 and 
:b=i6.58, and the green sphere (inside) with a=10.0 and b=63 js. lie images are rendered usmg 
Ray-Citing. 
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(b).exfemal sphere 

V Figure 8. The segmettfatioinof control 3D Wdboll distributed volume data. The color is 
. " .E©B(155,^0)/w1bere a' is file value at the point (r, 
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(a) 




(b) 



Figure 9- Test bed foi our Weibull E-SD modeling scheme. They come trom two different 
€3q)eriinents. The siZiC of image (a) is 5 12*5 1 2* 1 24, the grayrlevel is from Q to 255, and 
there is only one cell, its spmdle is at the up part of the cell. The size of image (b) is 
5 12*5 12* 146, has the same gray-level with (a) 
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3.2 Laser Confocal Microscope Data 

Li this example we use data ffoin two different LCM ejcperiments designed for 
investigating the meiotic spindte wi&in a mouse celt, and segmenting, tlie spindle. 

Figure "9 shows a LCM test bed for our Weibitil E-SD modeling scheme from 
different experiments. The size of data in Figure 9(a) is 512*512*124, and the gray-level 
is from 0 to 255. The spindle is at the upper left and bri^ter part of the cell. The size of 
data in Figure 9(b) is 512*512*144 .and his flie same graj^-Ievel range as in Figure 9(a). 
In Figure 9(b) the egg pictured has a meiotic spiiidie at the Upper left and the remains of 
the primary polar body at the bottom right 

Figure 10 shows the E-SD plots of Figure 9. The colors have the same meanings 
as in Figure 7. Li Figure 10(a), the size of the window is (170, 237; 4, 27) corresponding 
to the segmentation FigufelKaXJh Figure rO(b)V the size of the window is (183, 255; 2, 
30) correspondmg to the se^eflfeficaj Figure n(b). We set the threshold 2;=34, and h=2- 
Figure II shows the spindle segmentation of FJgiire 9. llfesegmentation in Fi^^ II (a) 
inchides 3902 voxels, andai308 voxels ill Figii^ 

What gives nse to i^chMQar^Jmgers hii^'-f^^ not Jmown. A piausibie 

explanation is that a jfirt^er corraspotiia^.:^: att SDO. When we move the window on 
FigurelO (a) to the smallfinger, it segments tiie region ofpart of the membrane of the cell 
in Fig 9 (a) (see Figure 12), 
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(a) The E-SD plot of image (a) of Fig 9. 




(b) The E-SD plot of image (b).Gf Fig 9. . 



Figure 10. The E-SD plots of Fig 9, respectively. The colors have the same means as Fig 
7. In (a), the size of window is (178, 237; 4, 27) corresponding to the segmentation (a) 
in Fig 11, and the threshold T,=34, so fliere is. a blank on the left side. There are two dear 
fingers. In (b), the size of window is (li67, 255; 2, 30) correspdhding to the segmentatJob 
(b) in Fig 11, and the threshold r,=34, so there is a blank on the left side. ThetB are two 
clear fingers. The box size is 2*2*2 




Figure 1 1. The spindle segmentation of Fig 9. respectively. The color is RGB(0, v,0}, whei« 
v is the value at. the point<x.y,r). ft arff i^idetedby IK 

segmentation in (a) indgd«( 3902 boxes,, and 11308 boxes in (b). 
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(a) (b) 
Figure 12. The region corresponding to the small finger in Fig 12 (a). It is rendered by 
using OpenGL point clouds from two different views. 



4. Conclusion 

We have proposed a coarse-grain approach for the segmentation of an object from 
volume data based on the WeibuII E-SD field. We have shown that ohe can decide tlte 
noise index in a ^-voxel by using the WeibuII law, and use the E-SD plot as a guide to . 
segment an object. We have consistently demonstrated this approach on a controlled as 
well as real volume data. 

Our future work includes finding out the relationship between volume data and its E- 
SD field, and applying our model to other types of volume data. 
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WHAT IS CLAIMED IS: 

1, A computer system for storing, archiving, querying and retrieving information 
relating to tliree-dimensional objects; the system comprising: 

data acquisition means for acquiring point coordinate data about the three-dimensional 

object; 

a database.component; 

a processor operable to: 

generate modeled data from the point coordinate data; , ^ 

segment the modeled data into feature daita represaitirig a plvurality of features of 

the object; 

store the modeled data and the feature data in the database Component; and 
retrieve modeled data and feature data from the database using search criteria 

conlprising representing an object feature; and 

a user interfece operative with flie processor to: 

input to tiie processor search criteria; and ' 

display data retrieved by the processor as a representation of an object feature • 

2. The system of claim 1 wherein the point coordinate data is surface daita. 

3 . The system of claim 1 wherein the point coordiaate data is volimie data. 

4. The system of claim 1 wherein the feature data represents a point. 

5. The system of claim 1 wherein the feature data represaits a curve. 

6. The system of claim Iwhdrein the feature data represaats a facet on a surface. 

7. The system of claim 1 whaein the processor is furflier operable to compress the 
modeled data, 

8. The system of claim 7 wherein the modded data conq)rises a tiiangie mesh and 
■:i the prbcessof is operable to corapress fee modeled data using B-spi^^ 

9. The^stiktiofcteifflTwhei&infeemo^^^ 

. ;i fee processdr is operable fo segnefflt fe^ inodelisd data usmjg subdiviaoii siafece canpressioiL 

10. ■ -The- systoii of claiia I wterein '-fee poiiit coordinate dabi cpinpii^es a triangle , 
mesh and the processor is operable to segment fee mod^ed xlata using watrashed segmentatiop 
mefeod including improved curvature estimation. 

11.. The system of claim 1 wherein fee poiirt cwordin^ 
v ffltesh^ fee jM:dcesSor is oper^le to segment fee niodeled data nang a hybrid segmentation 
mefeod. 

: :'^^-'t^/-, 12. The sj^tem^^rf claim 1 v;^^^ pdnt cootdioate data compdsos a triangte 
'rxi;^,f mesh and fee processor is opoable to segment fee ikodeled data uang watershed segmentation. ^ 
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13. The system of claim 1 wherein the point coordinate data comprises volume data 
the processor is operable to segment the modeled data using Weibull E-SD fields. 

14. The system of claim 1 wherein the point coordinate data comprises volume data 
and the processor is operable to segment the modeled data using Greedy connected component 
labeling refinement. 

15. The system of claim 1 wherein tlie user interface includes a graphic user interface 
that operates with the processor to allow a sketch-based search of the database. 

16. A method for storing, archiving, querying and retrieving information relating to 
3D objects; the system comprising: ' 

acquiring point coordinate data from a physical object; 
generating modeled data from the point coordinate data; 

segmenting the modeled data into feature data representing a plurality of features of tible 
. physical object; 

storing the modeled data and the feature data in a database; and 
organizing the modeled data and the feature data so that features of the physical "object 
can hp automatically extracted for online query and retrieval of fiie plurality of features of the 
physical object 

17. The inethod of claim 16 wherein the point coordinate data is surface data. 

18. The method ofclaini 16 wherdn the point coordinate data.is volume data. 

19. The method of claim 16 further comprising compressing flie modeled data. 

20. The method of claim 19 wherdn Ae modeled data comprises a triangle mesh aad 
. cpinpressing (he modeled <kta is performed uang B-spline curves- 

21. The method of claim 19 wherein the modeled data comprises a triangle mesh and 
conipressiiig the modeled data-is performed using a subdivision surface compression method. 

22. The method of claim 16 wherein the point coordinate data comprises .a triangle 
, mesh and segmenting the modded data is performed using watershed segment^on, 

. . 23. The method of claim 16 wherein the point coordinate data conpises a triaQgle 
. mesh and. segmenting the modeled data is performed using a hybrid segmentation rtiethod. 
., , . 24. ■ The method ofelaim 16 wherein the point coordinate data comprises volume data 
-and segmenting the modeled data is performed using WeibuU E-SD fields; 

. 25. : The method ofehiim 16 wherdn the point coordinate data <^ 
aiid segmenfing the modeled .d^ is performed Bang Greedy amected componeht labeling 
refinanent 
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